{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Distance inference between contigs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By now we have inferred distance between contigs based on Maximum Liklihood Estimate. The results are in file `distance.txt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#Contig1\tContig2\tLength1\tLength2\tLDE1\tLDE2\tLDE\tObservedLinks\tExpectedLinksIfAdjacent\tMLEdistance\n",
      "jpcChr1.ctg199\tjpcChr1.ctg257\t124567\t274565\t0.3195\t2.0838\t1.1607\t2\t27.4\t1617125\n",
      "idcChr1.ctg353\tidcChr1.ctg382\t143105\t270892\t2.1577\t1.0544\t1.3505\t2\t34.2\t2190000\n",
      "jpcChr1.ctg373\tjpcChr1.ctg379\t153398\t63193\t1.2697\t0.4859\t0.9594\t2\t12.0\t325545\n",
      "jpcChr1.ctg125\tjpcChr1.ctg174\t178182\t134096\t1.6988\t1.2466\t1.4874\t2\t30.4\t1482910\n",
      "idcChr1.ctg385\tidcChr1.ctg461\t201792\t195332\t2.2693\t2.6181\t2.4346\t5\t63.6\t1548564\n",
      "idcChr1.ctg104\tidcChr1.ctg146\t215933\t64665\t1.1102\t0.8443\t1.0423\t2\t15.0\t524288\n",
      "jpcChr1.ctg67\tjpcChr1.ctg76\t99865\t11957\t1.6615\t0.2000\t1.3249\t2\t4.1\t37315\n",
      "idcChr1.ctg284\tidcChr1.ctg48\t218292\t301915\t0.6696\t2.0850\t1.2945\t2\t43.1\t3683127\n",
      "idcChr1.ctg314\tjpcChr1.ctg305\t50575\t282724\t0.2000\t1.7503\t1.2594\t3\t16.8\t387141\n"
     ]
    }
   ],
   "source": [
    "!head distance.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The truth is in file `ctg.posi.bed`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "idcChr1\t1\t127964\tidcChr1.ctg1\n",
      "idcChr1\t127965\t176555\tidcChr1.ctg2\n",
      "idcChr1\t176556\t296400\tidcChr1.ctg3\n",
      "idcChr1\t296401\t481842\tidcChr1.ctg4\n",
      "idcChr1\t481843\t725279\tidcChr1.ctg5\n",
      "idcChr1\t725280\t726201\tidcChr1.ctg6\n",
      "idcChr1\t726202\t811858\tidcChr1.ctg7\n",
      "idcChr1\t811859\t1016514\tidcChr1.ctg8\n",
      "idcChr1\t1016515\t1090313\tidcChr1.ctg9\n",
      "idcChr1\t1090314\t1122899\tidcChr1.ctg10\n"
     ]
    }
   ],
   "source": [
    "!head ctg.posi.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "matplotlib.rcParams['ps.useafm'] = True\n",
    "matplotlib.rcParams['pdf.use14corefonts'] = True\n",
    "plt.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Helvetica']})\n",
    "#matplotlib.rcParams['text.usetex'] = True\n",
    "\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>#Contig1</th>\n",
       "      <th>Contig2</th>\n",
       "      <th>Length1</th>\n",
       "      <th>Length2</th>\n",
       "      <th>LDE1</th>\n",
       "      <th>LDE2</th>\n",
       "      <th>LDE</th>\n",
       "      <th>ObservedLinks</th>\n",
       "      <th>ExpectedLinksIfAdjacent</th>\n",
       "      <th>MLEdistance</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>jpcChr1.ctg199</td>\n",
       "      <td>jpcChr1.ctg257</td>\n",
       "      <td>124567</td>\n",
       "      <td>274565</td>\n",
       "      <td>0.3195</td>\n",
       "      <td>2.0838</td>\n",
       "      <td>1.1607</td>\n",
       "      <td>2</td>\n",
       "      <td>27.4</td>\n",
       "      <td>1617125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>idcChr1.ctg353</td>\n",
       "      <td>idcChr1.ctg382</td>\n",
       "      <td>143105</td>\n",
       "      <td>270892</td>\n",
       "      <td>2.1577</td>\n",
       "      <td>1.0544</td>\n",
       "      <td>1.3505</td>\n",
       "      <td>2</td>\n",
       "      <td>34.2</td>\n",
       "      <td>2190000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>jpcChr1.ctg373</td>\n",
       "      <td>jpcChr1.ctg379</td>\n",
       "      <td>153398</td>\n",
       "      <td>63193</td>\n",
       "      <td>1.2697</td>\n",
       "      <td>0.4859</td>\n",
       "      <td>0.9594</td>\n",
       "      <td>2</td>\n",
       "      <td>12.0</td>\n",
       "      <td>325545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>jpcChr1.ctg125</td>\n",
       "      <td>jpcChr1.ctg174</td>\n",
       "      <td>178182</td>\n",
       "      <td>134096</td>\n",
       "      <td>1.6988</td>\n",
       "      <td>1.2466</td>\n",
       "      <td>1.4874</td>\n",
       "      <td>2</td>\n",
       "      <td>30.4</td>\n",
       "      <td>1482910</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>idcChr1.ctg385</td>\n",
       "      <td>idcChr1.ctg461</td>\n",
       "      <td>201792</td>\n",
       "      <td>195332</td>\n",
       "      <td>2.2693</td>\n",
       "      <td>2.6181</td>\n",
       "      <td>2.4346</td>\n",
       "      <td>5</td>\n",
       "      <td>63.6</td>\n",
       "      <td>1548564</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         #Contig1         Contig2  Length1  Length2    LDE1    LDE2     LDE  \\\n",
       "0  jpcChr1.ctg199  jpcChr1.ctg257   124567   274565  0.3195  2.0838  1.1607   \n",
       "1  idcChr1.ctg353  idcChr1.ctg382   143105   270892  2.1577  1.0544  1.3505   \n",
       "2  jpcChr1.ctg373  jpcChr1.ctg379   153398    63193  1.2697  0.4859  0.9594   \n",
       "3  jpcChr1.ctg125  jpcChr1.ctg174   178182   134096  1.6988  1.2466  1.4874   \n",
       "4  idcChr1.ctg385  idcChr1.ctg461   201792   195332  2.2693  2.6181  2.4346   \n",
       "\n",
       "   ObservedLinks  ExpectedLinksIfAdjacent  MLEdistance  \n",
       "0              2                     27.4      1617125  \n",
       "1              2                     34.2      2190000  \n",
       "2              2                     12.0       325545  \n",
       "3              2                     30.4      1482910  \n",
       "4              5                     63.6      1548564  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(\"distance.txt\", sep=\"\\t\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[0;33m14:22:14 [base]\u001b[0m\u001b[0;35m Load file `ctg.posi.bed`\u001b[0m\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "idcChr1\t1\t127964\tidcChr1.ctg1\n",
      "idcChr1\t127965\t176555\tidcChr1.ctg2\n",
      "idcChr1\t176556\t296400\tidcChr1.ctg3\n",
      "idcChr1\t296401\t481842\tidcChr1.ctg4\n",
      "idcChr1\t481843\t725279\tidcChr1.ctg5\n",
      "idcChr1\t725280\t726201\tidcChr1.ctg6\n",
      "idcChr1\t726202\t811858\tidcChr1.ctg7\n",
      "idcChr1\t811859\t1016514\tidcChr1.ctg8\n",
      "idcChr1\t1016515\t1090313\tidcChr1.ctg9\n",
      "idcChr1\t1090314\t1122899\tidcChr1.ctg10\n"
     ]
    }
   ],
   "source": [
    "from jcvi.formats.bed import Bed\n",
    "\n",
    "bed = Bed(\"ctg.posi.bed\")\n",
    "print \"\\n\".join(str(x) for x in bed[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = []\n",
    "order = bed.order\n",
    "for i, row in df.iterrows():\n",
    "    at, bt = row[\"#Contig1\"], row[\"Contig2\"]\n",
    "    pair = \";\".join((at, bt))\n",
    "    nlinks = row[\"ObservedLinks\"]\n",
    "    if nlinks < 2:\n",
    "        continue\n",
    "    inferred_dist = row[\"MLEdistance\"]\n",
    "    ai, abed = order[at] \n",
    "    bi, bbed = order[bt]\n",
    "    if abed.seqid == bbed.seqid:\n",
    "        if abed.start > bbed.start:\n",
    "            abed, bbed = bbed, abed\n",
    "        true_dist = bbed.end - abed.start\n",
    "    else:\n",
    "        true_dist = 100000000\n",
    "    data.append((pair, nlinks, inferred_dist, true_dist))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Pair</th>\n",
       "      <th>ObservedLinks</th>\n",
       "      <th>InferredDist</th>\n",
       "      <th>TrueDist</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>jpcChr1.ctg199;jpcChr1.ctg257</td>\n",
       "      <td>2</td>\n",
       "      <td>1617125</td>\n",
       "      <td>4852517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>idcChr1.ctg353;idcChr1.ctg382</td>\n",
       "      <td>2</td>\n",
       "      <td>2190000</td>\n",
       "      <td>3481690</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>jpcChr1.ctg373;jpcChr1.ctg379</td>\n",
       "      <td>2</td>\n",
       "      <td>325545</td>\n",
       "      <td>632562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>jpcChr1.ctg125;jpcChr1.ctg174</td>\n",
       "      <td>2</td>\n",
       "      <td>1482910</td>\n",
       "      <td>5281491</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>idcChr1.ctg385;idcChr1.ctg461</td>\n",
       "      <td>5</td>\n",
       "      <td>1548564</td>\n",
       "      <td>6924766</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            Pair  ObservedLinks  InferredDist  TrueDist\n",
       "0  jpcChr1.ctg199;jpcChr1.ctg257              2       1617125   4852517\n",
       "1  idcChr1.ctg353;idcChr1.ctg382              2       2190000   3481690\n",
       "2  jpcChr1.ctg373;jpcChr1.ctg379              2        325545    632562\n",
       "3  jpcChr1.ctg125;jpcChr1.ctg174              2       1482910   5281491\n",
       "4  idcChr1.ctg385;idcChr1.ctg461              5       1548564   6924766"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf = pd.DataFrame(data, columns=[\"Pair\", \"ObservedLinks\", \"InferredDist\", \"TrueDist\"])\n",
    "tf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7334, 4)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.JointGrid at 0x1a1b2cc7d0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAGoCAYAAADmTPpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8HMXdBvBnTiederPVJVe5S7bljgsdY6rpLZQkpoSQkIRQk7z0JJAQkhAgNIPBMQ4QejOYZhsXjMFN7rItq1hW7/3u5v1DkpFlna5tvXu+n48wOt3tzp6kffSbnZ0RUkoQEREZgUXvBhAREfVgKBERkWEwlIiIyDAYSkREZBgMJSIiMgyGEhERGQZDiYiIDIOhREREhsFQIiIiw7Dq3QAPcdoJIgoEQu8GGJ1ZQomIemlo68SR+jbUtXSitdOB2HAr4iPDkJkQgdAQdoCQeTGUiEzA4ZRYW1CFj7aXYdOhWuyvaOq3+yA0RGBsaiymDk3AaeOSMXP4IIRZGVJkHsIkE7KaopFESmts68RLawux7JtDKG9oR1RYCEalxGBUcjTS4sIRHR6KsBALmjvsaGyzo6S2BQermlFQ0YR2uxPRNivOm5SOK6ZnYWJmHIRg75HO+A1wg6FEZEDtdgdeWluIf3+1H/WtnZiUFYeTRydjypAEjyqfDrsT20vr8c3BanxzsAYddidyMmJx04kjcVZOKqzs4tMLQ8kNhhKRwWwqrMFdb27D/spmTM6Kx6VTMzEiKdrn7bV02PF1QRVW5B9BWX0bshIjcNsZo7FwUgYsFp4jNcY33A2GEpFBtNsd+MuKPVj89UEkxdjw0znDMDkrQbHtO50S3x2qxVubS1BY3YKxqTG466yxOHl0Erv1tMM32g2GEpEBHKhswi+Xb8aOww2YPz4FV84YgvDQEFX25ZQS6/dX443vilHe0I6ZwxNxz9njMDkrXpX90TEYSm4wlIh0tiL/CG57fQtCLAI3nTgSU4cqVx0NxO5w4ovdFXhrcynqWztxxfQs3LVgLBKiwjTZf5BiKLnBUCLSicMp8Y/P9uJfXxQgOykKvz59NAZF2zRvR2uHA29tLsFH28sQFxGKe84eh0unZrJLTx18U91gKBHpoL61E7/+72Z8uacSJ49Owk/mDNf9fqJD1c14ce1B7C1vwrShCfjTRbkYnRKja5sCEEPJDYYSkcb2ljfihlc2oaS2FdedMBSnj0sxTFXilBKr9lTi1Y1FaLc7cOeZY7Fo7nCO0lMO30g3GEpEGnp3SynufnM7bFYLfnX6KIxNjdW7Sf2qb+3EC2sOYNOhWswanoi/XT4ZGfERejcrEDCU3GAoEWmg3e7Awx/swtINhzA2NQa/PHUUEg0+oEBKia/2VuKV9YUIDbHgoYU5WDg53TBVnUnxzXODoUSkspLaFtyy7HtsLanHOblpuGJGFqwW88yoUN7Qhn9/tR97yhtxweR0/PmiiYgIU2e4ehBgKLnBUCJS0ee7ynHb61vR6XDiphNHYsbwRL2b5BOnU+KdLaX433clGJsWg+eumYasxEi9m2VGDCU3GEpEKqhv7cSD7+/Em9+XYOigSPzqtFFIizP/NZktxbV48ssCWC0WPHlVHuaNStK7SWbDUHKDoUSkICkl3t1yGH/8aBeqm9qxcHIGLsrLCKgJUI/Ut+Hxz/agtLYVdy4Yi5tOHMHrTJ7jG+UGQ4lIIZsKa/Dnj3bju6JajEyKwk/mDMdIPyZSNbK2TgeeXb0fGw7U4KIpGXj04olcXNAzDCU3GEpEfpBSYm1BNZ76sgDrD1QjLiIUl0/Pwkmjk2AJ8OpBSom3NnddZzppdBKe/tEURNm4bqgbgf1DoQCGEpEPDte14qPtZVj2TREOVjUjPjIU501Mx2njkmGzBtfItM93l+PFrw8iJz0OL/1kui5TJZkIQ8kNhhKRB6SU2HG4ASt3luOzXeXYcbgBADA6JRqnjU3BrBHBvez4pkM1+NfnBUiPD8fSRTM5Ms81hpIbDCUiF9rtDmw4UIPPdpZj5a5yHKlvgwAwJjUGeUMSMG1oAtI5y8FRe4404rFP9yAiNAT/uX4mxqRy3rx+MJTcYCgR9dLcbsdnu8rx6Y5yfLW3As3tDtisFkzMjMPUoQnIy0pAbESo3s00rJLaFvz5492QUuLVG2ZhXJoxp1HSEUPJDYYSBb0OuxOr91bi3S2lWLmzHG12JxIiQzFlSAKmDE1ATnpcUHfNeetIfRse/nAn7E6JZdfPRE5GnN5NMhKGkhsMJQpahVXNWLKuEG93L3IXE27FzOGDMGfkIIxOjQn40XNqKm/oCqYOhxPLFs1CbiaDqRt/qNxgKFFQkVLim4M1WPz1QXy2sxwhFoEZwxMxN3swcjPjTDUnndFVNrbh4Q93obXTgaWLZnK59S4MJTcYShQUnE6Jj/OP4OmvCrDjcANiwq04fVwKzhifgoRIY8/WbWaVje14+MOdaOlwYOmiGcgbos1S7wbGUHKDoUQBTUqJlTvL8fjKvdh9pBHp8eE4OycN80Yl8TqRRqqb2vHQhzvR3O7Af64P+oqJoeQGQ4kCUs9aQI9/uhfbS+uRGhuOi6dmYvaIQVxFVQc9wdRVMQV1MPGHzw2GUj/K6lux8WANdpY1oKapA3OyB+PkMUmIZzeP4UkpsW5/NR77dA82F9UhOcaGC/MyMG9UEkIYRrqqavqhK+8/i2ZiUnAGE38I3WAo9dLpcOL5NQfwj5X70OFwwmoRiAgNQWO7HSFC4JoThuLus8YiPDS4ppExi+8O1eCvK/Zgw8EaDIoKwwV5GTh5dFJAzdBtdlVN7Xjog51o7QzaYGIoucFQ6lZa14obX9mEHYcbMH1YAi6ekomMhAhYhMCByiZ8tacSn++uwLi0GDx51ZSAnf3ZjHYcrsffPtmDL/ZUIj4iFAsnp+PUsSm8ZmRQvSumZdfPxMTMoAomhpIbDCV0Ta552bPrUdvcgRsHWB30+6JaPLNqPwBg+Q2zeFOgzvZXNuHxlXvx4bYyRNusOHdiGs6ckMpK1gR6RuW1dXYNfgiiYGIouRH0oVRW3xVINU0duOfscW4roMrGdjz4wQ44nBJv/OwEZCdzfi+tldS24J+f7cOb35cgzGrBWTlpOCc3jcsmmEzvYFp2fdDcYMtQciOoQ6mlw44Ln1qH4toW3HPWWI8Dpqy+FQ++vxOhVgveunk2Z0TWSEVjG57+cj/+s+EQhADOGJeC8ydnII5z0ZlWEAYTQ8mNoA0lKSV+uXwzPtpehrsWjPW6+6C4pgUPfLAD6XEReOvnsxETzhOjWupaOvDs6gNYsrYQ7XYHThmTjAvzMrhuT4CobGzDQx/uQofdif8smhnowcRQciNoQ+n51Qfwx4924YrpWVg4OcOnbWwvrccjH+/CSaOT8MJ10znkWGFN7Xa89PVBPLv6AJrb7ZidPRiXTMlEaly43k0jhfUEU2uHA4uvm4aZIwbp3SS18CThRlCG0rqCKly9+BtMH5aIX502CsKPiTdX7izHi2sPYtHc4fi/c8cr2Mrg1dbpwH82HMLTX+1HTXMHpg9LwCVTszCE3aQBrbqpHX/+eDeqmtrx1FVTcPr4FL2bpAaGkhtBF0oltS04919fI9pmxYPn5yAizP+RWkvWFeKTHUfw54tyceWMIQq0Mjg1tduxbMMhPLf6AKqbO5CbEYfLpmUhO5nD74NFQ1sn/rJiNwqrWvDAwgm4etZQvZukNIaSG0EVSm2dDlz873U4WNWMhxfmIE2hVUMdTom/frIbOw43YOmimThhZMB2PaiivrUTL68rxOKvD6K+tRMTM+NwweQMLhAXpNo6HXjii33YXFSHRXOH43dnjwukrvGAORC1BE0oSSlx+xvb8Ob3Jbhj/hhMGarsbMUtHXbc994ONLXZ8c4tczBscJSi2w9Eu8oasOybQ3jr+1K0dDgwdUgCLsjLYGVEcDollm44hBU7juDk0Un4xxWTA2WaL4aSG0ETSq+sL8S97+7AxVMycMnULAWadLzyhjbc+24+Bkfb8PYtczhUuR9tnQ58tL0My74pwneHahEaIjBrxCCcnZuGYYMY5HSslTvL8cr6QqTGhuOZa6YGwg3rDCU3giKUNh6swVXPb8DEzDj8dv4YVVcU3VXWgD99tAuzRgzCkp9M57xrAFo7HFhbUIWP8suwcmc5GtvsSI8Lx2njUnDiqCREh/OmV3KtoKIR//hsHxraOnHHmWOwaO4IM3fnmbbhWgn4UDpQ2YRLnlkPm9WChy/IQWSY+ifAL/dU4LnVB3BRXgb+csnEoAumysZ25JfWY3tpPb45UI2NhTXodEhE2UIwbWjXKq8T0mP9GvVIwaWhrRMvrDmAbwtrMWNYIv566UQMNWdlzR96NwI6lMrqW3Hx0+vQ1GHHfedOQLpCAxs88fbmUry+qRgLclLxzysmw2YNrPnYpJSobu5AYVUzCqtbUFjVjL3ljdhWUo8jDW0Aun77shIjkZsRh4mZcRifFht0AU3KkVJi9b5KvLzuEBxOiVtOycZNJ40w21yHDCU3AjaUKhvbccVz61FW34Y/nDMew3UYePBxfhleWX8Ic7MH44kr85AYZb4LtR12Jw5VN2N/ZRP2VzZjf0UTCiqbcKCyGU3t9qPPswggNS4cwwdFYfjgaAxPisKwQZGaVKYUXGqaO/CfDYew/kA1MhMi8Nv5o3H+pAyzdOmZopF6CshQ2lfeiJ+89C0qmtpx94Kxug4tXrW3Ai+sOYjEqDD844rJmD1ysG5tGUhdS0dX8FR0BVBBd/iU1LTC0etnZFBUGNLiI5AeF460uHCkxIYjNS4cSTE2WC2sgkg720vrsXxjEQ5WNWN0SjRuPnkkzp2YjlBjV+MMJTcCLpRW763ELa9+jxAhcPuZYwyx7tHBqmY8+cU+lNW34dJpmbj1tFHITNB+doIOuxMltS04VN1yXOVT09xx9HlWi0BaXDjS4yN++IgLR1pchCI3GxMpxSklvjlQgzc3l6C0thUpsTZcM2soLp6aibQ47brrvcBQciNgQqm+pRN/+mgXXttUjKyECNxx5lgkxRhnws62Tgde31SMlTvLAQAX5mXg/MnpOGHEIMWus3Q6nKhqakd5QzsO17XiUHULimqacai6K4jK6lvh7PVOxoRbuwMnAunxXSGUER+BpGgbLOboCiEC0BVOW4vr8FF+GfJLG2ARwNzswThnYhpOH5dipMl7+YvlhulDqb6lE69uLMLirw+gprkD5+Sm4ZKpWYZddbS6qR3vbCnF2oIqtHY6ERcRiomZccjJiMOwQZEYHG1DfGTX/U1S/nDgHXYn6ls7j340tHaiuqkD5Y1tqGhoR3lDG2qaO457o2LDrUiJDUdyjK3r39hwpMTakB4fgVjObE4BqLyhDav2VmJtQRUqGtthEUBORhxmDk/EzOGDMH1YIuIidfvZZyi5YcpQqm3uwJqCKny5uwIrdhxBa4cDORmxuHL6EIwwQHedJzrsTmwrqcN3h2pRWN2M4tpWOJyefy9CLAKx4VYkRIUhPiIMCZGhiI8MQ0JUKBIiwpAQFYaUWBsHGlDQklLiUE0Lvi2swc7DDSioaILdKSEAjEmNwbi0WGQnR2NkUjSyk6MwdFCUFtejGEpumCKUnvxinyxvaEdJbQt2H2lEWX3XkOOYcCumDEnAgpxU088G0Olwoq6lqwpq7h7V1nMbjxACVotAlM2KqLAQRNmssFktvM+HyAsddicKKpuwu6wBe8sbUVLbiupe11ItAhgUZUNSzA8fg6LCun7vev3uhYeGwGoRCLEIWEMErBZL1/9bBCxCwCElHE4Jp5SwO7r+dTglHFLilDHJ/KV1wxShNOzuD2VMuBXJMTYMGxyF4YOjMDEzDmNTY80yDJSIDKilw47imlYU1bSgpLYFNc0dqGnuQF1L59F/OxxOxfZX+Mg5PGG5YYpQEkKsAGDMsdTHGgygSu9G6CRYj53HHVz8Pe4qKeUCpRoTiEwRSmYhhNgkpZymdzv0EKzHzuMOLsF63Foy5hA1IiIKSgwlIiIyDIaSsp7TuwE6CtZj53EHl2A9bs3wmhIRERkGKyUiIjIMhhIRERkGQ4mIiAyDoURERIZhitk6Tz3jTPn6ux/q3QwiIr8MjrR6NM1QxIipcsWKFThpdJLaTdKSR8duikqppjoYZzMhomDlbG1AaW2r3s3QhSlCiYgomAgARTUtejdDFwwlIiKDCQ2xYFdZg97N0AVDiYjIYGyhFmwrqUMwTm5gioEOFBwc9k7UVZShs6O9ay14IjMSAqFhNsQnpyHE6tuy6+HWENS2dGJveRPGpMYo3EBjYyiRYdRVlCE+NhYJiYlcVZdMS0qJ2ppq1FWUYVD6EJ+2EWkLQQeAT3ccCbpQYvcdGUZnRzsDiUxPCIGExEFdFb+PrBYLRiVH4/1th4OuC4+hRMYhJQOJAoIQwu8u6FPGJmNveRNW7wuuW2IYSkREBjQ3ezASo8Lw7y8LgqpaYigRBalPVqzAhPFjMW7MKPzl0Uf6fc4rLy9Bemoypk3Nw7SpeXhx8QsAgC1btmDenNmYNDEHU/Im4fXXX9Oy6UetWb0aM6ZPRYQtFG+++T+vX//qq8swJW8SpuRNwolz52Dr1q0AgOLiYpxx2qnIzRmPSRNz8K8n/nncax//22MIs1pQVaVOJRMaYsF5E9Ox4WAN3t9Wpso+jIgDHYhUZrfbYbUq/6vmcDgQEhLi82t/desv8NGKT5GZmYkTZs3Aueedj/Hjxx/33Esvuwz/fOLJYx6LjIzEi0texqhRo3D48GHMmjEN8+efifj4eJ/a46usIUPwwuKX8PfH/+bT64cPG47Pv/gKCQkJWPHxx/j5z27C2vUbYLVa8Ze/Poa8KVPQ2NiImTOm4bTTzzj6/hQXF+Pzzz7DkCG+DWTw1PzxKVi7vwr3vpuP2SMHYXC0TdX9GQErJaJuhYWFyJkwDj/9yY8xJW8SLr/sUrS0dN1V//133+G0U07GzBnTcM5ZC1BW1vWX6+IXnscJs2Zg6pTJuOzSS44+f9FPf4I7fnsbzjjtVPzu7ruwetWqo9XG9GldJzopJe6+8w5MnpSLvMkTj1Ybq776Cqefegouv+xS5EwYh2uvufpo982okcPx8EMP4uQT5+F//3vD52P9duNGjByZjREjRiAsLAyXXXY53n/vXY9fP3r0aIwaNQoAkJ6ejqTkZFRWVgIA7r/vXrz//nvHvebBB+7Hj6+7FvNPPw3jx47G4hee97n9PYYNG4aJEyfCYjn+VPa3x/6KE2bNwJS8SXjg/vv6ff0Js2cjISEBADBz1iyUlpYAANLS0pA3ZQoAICYmBmPHjsPh0tKjr7v9t7fhT488qvo1UItF4KYTR6CpzY5f/XczOh1OVfdnBKyUiHrZu2cPnnvuBcyeMwc3XP9TPPPvp/HLW3+FX//qVrz59jtISkrC66+/hnv/7/d4/oUXccGFF2HR9TcAAO79vz/gpRcX45Zf/BIAsG/fPqz4dCVCQkJwwcLz8cQTT2L2nDloampCeHg43n77LWzduhXffb8FVVVVmD1rBubNOxEAsGXLZmzZlo/09HScNG8u1q1dizlz5wIAwsPD8dXqNce1/dVXl+Hxvz123OMjR2bjtdePDbDSw6XIzMo8+nlGZia+3fhNv+/J22+9hTVr1mDUqNF47G+PIysr65ivf7txIzo6OjBy5EgAwP0PPOjy/d2+fRu+Xrsezc3NmDFtCs46+xykp6cf85xTTjoRjU2Nx7320Uf/itNOP93ltntb+emnKCjYh3Xrv4GUEhddsBBrVq/GvBNPdPmal15cjDMXLDju8cLCQmzdshkzZs4EALz//nvIyEjHpEmTPGqLvzITIrFo7nA8u/oA7n5zOx67dGJADwhiKBH1kpWVhdlz5gAArrrqajz55L8w/8wF2LEjH2ctmA+gq+srLTUNALAjPx/33ft/qKuvQ3NTE86YP//oti665JKj3WuzZ8/GHbf/FldedRUuuPAiZGZmYt3XX+PyK65ASEgIUlJSMO/Ek7Bp07eIjYnF9OkzkJnZFRqTJk9C4aHCo6F06WWX99v2q676Ea666kceHWd/F877O9Gdc+55uPyKK2Gz2fDcs89g0U9+jE8/+/zo18vKyvDjH1+LF19c0m+10td5552PiIgIRERE4KSTT8G3327EwoUXHPOcL1et9ugYBvLZyk/x2cqVmD6tq9ppbmpCQcE+l6H01Zdf4qWXXsRXq44N+6amJlx+2SV47PG/IzY2Fi0tLXjkT3/CRys+8buN3jh5TDKqmtrx5vclSIqx4a4FYwI2mBhKRL30/UUXQkBKifHjJ2DN2nXHPf/6RT/BG2++jUmTJuGVl5dg1apVR78WFRV19P/vvOtunHX2OVjx8UeYN+cEfPzJSki4HlFls/1w7SAkJAQOu/247bbaj+3KeX35q3ji78dfWxk+ciSWLn/9mMeSUtNxqKj46DYKi4oxOCXtuNcOGjTo6P8vuv4G/O6eu49+3tDQgIXnn4sHHnwIM2fNcnksvfX3/valRKUkpcSdd92NG2686ZjH//30U1jcPVjjvfc/RHp6OrZt24af3XQD3vvgo2OOt7OzE5dfegmuvPIqXHjhRQCA/fv3o7DwIKZNmQwAKCkpwczpU7F2/TdITU31qG2+unhKJupaOvHMqv2oa+nAwxfkwBoSeFdgGEpEvRQVFWHD+vWYdcIJeO215ZgzZw7GjBmDqqrKo493dnZi7969mDBhAhobG5GWlobOzk4sf/VVpGdk9Lvd/fv3Izc3F7m5udiwYQP27NmNufNOxAvPPYdrrr0ONTU1+HrNajzy6F+wbcdOOKQ8Ghh2J9Dh6PpcSqDN7jwukADgsiuvwmVXXuXRcU6ZNh37CwpQePAg0jMy8NYbr+OFl5cet90jZWVITesKq/fffRejx45Fq92JEKcdl158Ea6++hpccsmlx7zm97+7B9NnzMAFF1x43H7ff/893HX3PWhubsbqVV/hj3/683HPUaJSOmP+mbj/vntx5VU/QnR0NEpLSxEaGoqbf34Lbv75LUefV1RUhMsvvRgvLXkFo0ePPvq4lBI33nA9xo4bi1//5rajj+fm5qK0rPzo56NGDsf6b77F4MGD/W6zO0IILJo7HDHhofjvt8WoaGzHk1flITIssE7jgXU0RH4aO24cli59GT//+c+QnT0KN/3sZoSFhWH5a2/gtl//CvUN9bDb7bj11l9hwoQJuP+BBzF39iwMGTIUObk5aGxs6ne7/3riH/jyy68QEhKCMWPH4aTTz0RYWBi+XrcOU/ImQwiB+//4COIGJwPYqfpxWq1WPPaPf+Ki886Gw+HA1df9GOPGTwAA/PGB+5A3dRrOPvc8PPPUv/Dxhx/AarUiISEB/37+RQDAq/99DWvWrEZVdTVefvllAMDTzy/GzKlTsCM/H+edd36/+50+fToWnncuiouL8Lvf/+G460ne2vTtt7j0kotQW1uLDz94Hw8+cD+2bsvHGfPnY/fuXZg3dzYAIDoqGkteWYrk5ORjXv/Hhx9EdXU1fvnLW46+Lxu++Rbr1q7Fsv8sRU5uLqZNzQMAPPTQH3HW2Wf71V5/CSFw+fQsJESFYsnaQlzw1Fo8/aMpyE4OnKmIhBluypo8Zar87Ov+L8JS4DhycC/GjB2n2/4LCwtxwcLzsGXrdp9e31/1EowuPPcsvP3Bx0c/j7B2dTE9+MD9iI6Oxm2/vV2vpmlqz+5dSB0++pjHPF15dvzEPLn47c8GfM7W4jo8vaoAnXaJhy7IwSVTMwd8vgEEzsqzREbR2t111t8HdekdSMAP75ndafw/gM1kUlY8/nzhRAwfHIXb39iKX/93M2qbO/Rult/YfUfUbdiwYVj/3daACpjWTomIUGOM0rrn/7ruFer9/vZUUeSbxKgw/P7scXhrcyne2VKK1fuqcN9543H+pHTTjs5jKJFxdI90U/uXyeyh09rpXcXh7fMBaBZkPd+LQAsnKSWgUShYLAKXTM3EjOGJeH7NAfzqv1vwzuZSPHxhLjLiIzRpg5IC6yeBTC00zIbammrFJp80ezdba6fs9yMQ9222781AetZTCg3TdkqgIYmReOC8Cbj2hKFYf6AaZzy+CkvWHoTDZN2mrJTIMOKT01BXUYbKqiqfpv3vNNkvX48Ohznb3SMsRPmKINRizq4nAMesPKs1i0XgrJw0TBuagMVfH8T97+/EO1sO49GLJ5pmsUCGEhlGiDXUq5U6t5f7voia1vLL2/RugmZyUsIV21ZuSuBPQKqGpJhw3LVgLNbur8Yr6wtxzhNr8POTR+KWU7Nhs/o2ia9WGEpkOkYPIz0CaEeFb+/JhGTlT/q9j9/fgOr5XjOcvCeEwNzswZiYEYelGw7hiS8K8MH2Mjx68URMH5aod/NcYiiRaRgxjMwUQN5sS6mw6nl/lAgnBpNvYiNCccsp2ZibPRiL1x7Epc+sx9WzhuCes8Yhyma8CDBei4j6MFoYaR1ESoaQr/v0N6SUCCcGk38mZcXjLxdPxBubirFsQxG+3leFf105BbmZcXo37RgMJTIso4RRMISQO0qFlL/hxGDyT3hoCK45YRimDUvEU18W4MKn1+LOBWNw/dwRsBhkcAmHhJMh6R1I+eVtRz+0sKOi/eiHGfjbVn/eV71/NgLBuLRYPHLRROQNicefPtqNm5d9h7ZOh97NAsBKiQxIr5NOoFRE+RUDH0dOsnKj43qOwZfKKb+8TdGReuSd6HArfnP6aHycfwRLNxzCT5d8i+evnab7dSaGEhmK1oEUCEHkLoQ8eb6/QeVrOPkaTOzGU4YQAmfnpiHKZsVzq/fjmsXfYOmimboGE7vvKGiZPZDyK9q8DiS1t+XLMQbTPVxGddLoJNx62ihsLqrDXz/Zo2tbGEpkGFpWSVqeCJW+VqRkGKmxba2ui/HakrJmDh+EM8an4OV1hdhaXKdbOxhKFFT0GLygJLXCSOn9eHvcrJaMoWsBwTA8umK3bm1gKBEREQAgMsyKvKx47D7SqFsbGEpERHRUcowNNc0daGq367J/hhKRStSYV45IbfWtnbDqeCMtQ4mCCu+LIXLN6ZSTf1NKAAAgAElEQVRYd6Aap45NRrROw8IZShR0tAwmpaulnORwRW9+HWg//mCVaE5f7KlAXUsnLszL0K0NDCUyDC1vhsxJCdcsnCYk21QJJzVoFXpkPPsrm/DyukLMzR6M+RNSdWsHQ4kMReu79LWumpQMJyUDRMlt+XKM7FbVV3VTO/75+T4kx9jwxJV5CNHxmhKnGSLDyU2xaXpjZO8Tohb3y/SctJW6h6lvmHhyj5Fa1ZCW3XacZkgZh+ta8eePd6Gt04kXb5iFxKgwXdvDUCJD0jqYemgZUL1P4EreZKtH95s/YcQqST8HKpvw6IrdCA2x4L83zkJOhv5rKzGUyLD0CqYegRBQauOABvPacKAaz6zaj8HRNvzn+pkYPjhK7yYBUDGUhBAvAjgXQIWUMqf7sckAngEQDsAO4OdSyo1qtYHMT+9g6sGA+oGi18VYJWnOKSXe2FSCd7aUYuqQBPz7milIjjHO90HNSmkJgCcBvNLrsb8AeEBK+bEQ4uzuz09WsQ0UAHquHRghnAD9AgrQL6TUqIj8CSReT/JNS4cdT31ZgO+L6nDF9Cw8sHACbNYQvZt1DNVCSUq5WggxrO/DAGK7/z8OwGG19k+Bx2jhBBx/Yg2EkNKiS44VkvbK6lrxt5V7Ud7QhocWTsDVs4ZCCGMsgd6b1teUfg3gEyHEY+gajj7b1ROFEDcCuBEAMrOGaNM6MgUjhlMPvUPK6JQIo0Ctknqf81IzMhXd9pbiOvzri32wWS1YumgmThg5SNHtK0nrULoZwG+klG8KIS4DsBjA6f09UUr5HIDnAGDylKlSuyaSWfQ+ORkxoADtQ8rIWB0NrPc5b/zEPMXOeav2VuC51QcwJiUGz183DZkJkUptWhVah9J1AH7V/f9vAHhB4/1TgOr717NZQgoI/KBiGOnng22HseybIszJHoRnr5mm23x23tC6hYcBnATgKwCnAtin8f4pSPTXxcOg0h4DST9vby7F65uKcXZuKv5++WTDDWhwRc0h4cvRNbJusBCiBMB9AG4A8E8hhBVAG7r7T4m0YOagMltIMYz09XVBFV7fVIwL8zLw2KWTdJ02yFtqjr670sWXpqq1TyJvubpobrSwMkM1pUUQBeogByXtPtKAZ1ftx8zhiXj04ommCiSAMzoQ9csM16jchYDaoaVHNbS9vJ3BNIBOhxPPrtqP9PgIPHvNVIRZzTfnNkOJyANmGOnXF7vQgs+K/CM40tCOl386A/GR+k6s6iuGEpGXzFBFUfBp63TgnS2lOHVMEk4anaR3c3xmvtqOyGByU2xHP0h9/COgf9tL69HS4cD180bo3RS/MJSIFMRgIr1sLqpDtM2K6cMT9W6KXxhKRApj1UR6OFTdjClDExAaYu7TurlbT2RgDCf1sAvveK2dDsRFhOrdDL8xlIiIAkC73YmIUPOf0jn6jkhlRlmoUG393Rel5rB03rN0rJhwK6qaOvRuht8YSkTkNU9vzO37PN47pZ6kaBuKa1r0bobfGEpE5BElZojo2YZS4cRq6Qfp8RHYUlyHtk4HwkPNMflqf8zfAUlEqsovb1N8yiI1thnsxqTEwO6U2FxUp3dT/MJQIlKZWa8naREcDCbljE6NgQDwbWGN3k3xC0OJiI6jZVj4uy+zhr7Som1WZCVGYuNBc4cSrykRqchsJ0xWLuY2JjUGawuqYHc4YTXpTbTmbDURKU7PQGIYKmNsagxaOhzYVdaod1N8xlAiIoZCgBibGgsA2Gji60oMJSKVmKXrjoEUOBKjwpAUY8P3RbV6N8VnDCWiIGW0YdlGaouZZSVEYF85u++IyEQYAIErPT4CB6ua4XBKvZviE4YSUZBhIAW2wdE2dDok6lrMOQ8eQ4koiDCQAl+0retOn7rWTp1b4huGElGQYCAFh55575rb7Tq3xDcMJSKVGGmi0EAOJCO9z0bg7L6WFGIROrfENwwlogBnlkDishbKsHeHktViztO7OVtNRB4xSyCRcupbuwY4DI4O07klvmEoEamIXUvq4vt7vMrGdkSEhiAxiqFERAbCKik4VTS2IyMhAkLwmhIR9UOPv+bNFki8nqScyqZ2ZCVE6N0MnzGUiAKM2QKJlFXV1I6sxEi9m+EzhhKRBnjtQ3l8T4/X3G5Hc7sDmayUiMgIWCUFt6qmrpnpMxNYKRGRG/zLntRW2dgVSlkMJSLSm1mrJA5yUE5FY0+lxO47IvIAqyVl8H3sX1VTOyLCQhAfGap3U3zGUCIKAGatkkhZTe12JESGmvYeJYChRGR6DCTq0dLhQFyEeaskgKFERBQwOh1OhIWY+7Ru7tYTBTlWSdSb1WJBp8Ocy6D3YCgRmVSgBFKgHIcRhIYItHY69G6GXxhKRCbEEzn1Z1C0DYfrWiGleaslhhKRyQR7IHE4uGspMTa0251H71cyI4YSkYa2l/t3sgj2QKKBDRnUNZPDluI6nVviO4YSkUkwkMid7KRo2KwWrCuo0rspPmMoEWnEnyop0APJm+Pzt9oMZNYQC8alxeLzXRVwOs15XcmqdwOIaGBaBtKOCt9O+BOSeZ3HKOZmD8aTXxZg/YFqzMkerHdzvMZQItKAr3/dqx1IvobQQNthQOlr+rBERNusWPbNIYYSER3PSIGkVAh5sg9vwym/vI0zhisgzGrBaeOS8d6Ww9hV1oBxabF6N8krvKZEZEBKBdKOivZjPrSk9f7oB+dOTEdEWAge+3SP3k3xGkOJSEW+VEn+BpJeITRQW0hb0TYrzpuUjs93VeDL3RV6N8crqoWSEOJFIUSFECK/z+O/FELsEULsEEL8Ra39E+lNy0AyUhD1x6jtCmTn5KYhMyEC97y1HY1tnXo3x2NqVkpLACzo/YAQ4hQACwFMlFJOAPCYivsn0o23gZRf3uZ1IBk9iHwR6EPftRQaYsFNJ45ARWMbHvpgp97N8ZhqoSSlXA2gps/DNwN4RErZ3v0cc9WVRB7wJZC8YdYgMmObzS47OQbnT8rA65tK8L/vSvRujke0Hn03GsA8IcQfAbQBuF1K+W1/TxRC3AjgRgDIzBqiXQuJ/KDmjZ1qnNTzK/yrTHKSvRstt6OinUPGXeh9zkvNyFRsu5dMzcS+ikb84e3tmJAea/jReFoPdLACSAAwC8AdAF4XLtbtlVI+J6WcJqWcNmiw+cbaU/BR6xqSUpVRfkXbcR9KbpP80/ucl5Co3DkvxCLwi1OyEREWgkVLvkVFo7G/V1qHUgmAt2SXjQCcAJg4ZHpqBpKvlA4gT/fnjrtj4nUl5cVHhuH2+WNQ3dyBRUs2obXDuGsuaR1K7wA4FQCEEKMBhAEw78yBRFAnkHytjoxQubBqMqYRSdH4xanZyC+txy+Wfw+7w6l3k/ql5pDw5QDWAxgjhCgRQiwC8CKAEd3DxP8L4Dpp5tWoKOipFUgeb0vjashT7trib3ckJ2X1zbShifjxnGH4fFcF7nxzmyEnbVVtoIOU8koXX7parX0SaUWtqYM8PVmrEUBlZWUuv5aWlub19vIr2rweCEHqmz8+FU1tdrzxXQliw0Nx33nj4eLSvi449x2Rl/QMJH/DaKDg8fZ1ngTVQMHEkXj6uTAvA83tdixZV4jw0BDctWCMYYKJoUTkITXXQ3J78d/HMPI1hLzZti9VFOlLCIGrZw1Fh0PimVX7ERYicNv8MXo3CwBDicgtta9fDBRI3oaRmiE00D4HCiZ24xmTEAI/mTMMDqcTT3xRAGuIBbeeNkrvZjGUiAaiRCANVCUpEUh6BJGrNrgKJ2+DictYaMMiBK6fNwJ2p8TjK/ciNMSCm08eqWubGEpE/VCqOlIzkJQMo7aiY+ZNRviQHJ+2465q6ovXlfRnEQI/O3EkHE6JR1fsRpjVgkVzh+vWHoYSUR9aBNKAr3MTSL6GUd/g8eW5voYVGZvFIvDzk7Nhd0o89MFORIWF4IoZ+kzvxlAi6kWrQHJVJQ0USN6GkTch5O02BwonV9USry0ZW4hF4JenZOMx+x7c89Z2RNqsOH9Suubt4CJ/RN0CIZDaivKPfqhJ7e2TPqwhFvzm9NEYmxaD217bgs92lmveBoYSBb3t5e2mDiStgqi//bpihMEX5BubNQS3zx+DoYMi8fNl32NtgbYzwTGUKKgZebqasrKyAU/uegRRf22gwBMZZsVdC8YiJc6G61/ehC3FdZrtm9eUKGgpHUhqVEmu+BIG7cXbvX6NLSvXo7b0d42pv2tLvK5kHjHhofjdWeNw33s7sGjJt3jnljnISoxUfb+slIgUoHQguaqQvKmO2ou3H/PhC09fq3XFlJvCYeRaiI8Mw50LxqLd7sSPX9qI+pZO1ffJUKKgpOQ1JKXX/xkokNzxN4TcbZeCT0Z8BH5zxmgcqm7Brf/dDIfKM4szlCjoaH0fkjdVkq+BpFVoMJiC0/i0WFw3exhW7a3EPz/bq+q+GEoUVAItkIxUwXDQQ2A7bWwyTh6dhCe+KMAXu9UbKs5QoqBh5EByxV0g6cGb/fo6NJzz3hlP1wSuwzE0MRJ3vLENVU3qjFxlKFFQ0HvqIHf6O3mrEUhtJTuP+yDyVJjVgltOyUZDWyfu/N82qLFwOIeEU8DTI5D87bZTKpA8CZ3ezwnPHO/RdtuLt3s0XJwCT1ZiJK6YPgRLNxzCe1sPY+HkDEW3z0qJyANqBZK3PBqe7UcVZIbqicPB9bdgQipGJEXhoQ92or5V2WHiDCUKaGqvh9SXJ0ua9+ZNleQukMwQKBQYLBaBRXOGo6a5A/9QeDQeQ4kClpECyd9uO08CSUlKbY9LpQeuEUnROHlMMpauP4TimhbFtstQooBk9EDyxkCBZLbqiFMMBZaLp2TCIgT+vlK5aomhRNQPNbvsAM+rJHeBpCYtwq6/VWc5HNw8EqPCcMb4FLyzpVSxasltKAkhjlt+sL/HiIzC3ypJyWHf/nTb6RlIgOcj8bTEQQ7Gc3ZuGoQQWPz1QUW250ml9I6HjxHpTutAUrPbzhUzdddR4EuMCsPskYPw2qZiNLb5PxLP5X1KQojRAMYBiBNCnN/rS7EAWF+T4ZghkPytkowSSP3dozTQEukU2OaPT8GafVX4YFsZrpzhX0faQDfPTgBwEYB4AJf2erwRwE1+7ZXIYJQMJFcCJZA81XfkXd9BDt5eT2LXnXGNTIpGZkIEXv+2WL1QklK+DeBtIcRcKeXXfu2FSGVariDrLpA87bYL5ECi4CKEwLxRSVi+sQjFNS1+LQboyTWlc4QQsUIIqxDiEyFEuRDiKp/3SKQwLbvtfA0kXycmBfQJpIEGObDrjvozc3giAOCTHUf82o4noXSWlLIBwLkAKgDkALjLr70SGYRegeTtaDuzcdd11x923ZlbSmw4hg2KxIp89UMptPvfswEsl1JWAlB36UEiD/lTJZkhkAK1SurvehKZ3+SseGwuqvNrFJ4ns4R/JITIB+AAcIsQYjAA7TrwiVwweiD5y9NAchVoes3i7UuVRIEhNzMe72w5jHX7q3HmhFSftuG2UpJS3gHgVABTpZSdANrQNSqPKOD5Msquhz9Vkr+B5Cu9qiR23QWG0cnRCA0R2FRY4/M2BrpP6SQp5are9ygJIXo/pdTnvRL5SYsqyZNAUqPbzpNA8iSMlFzzyNPtcALW4GYNsWDY4ChsLqrzfRsDfO0MAKtw7D1KPSSA93zeK5FO9Aqk/vha5ag1IMLbaYU8qZLYdRd8spOi8eXuCnQ6nAgN8X561YHuU/pD97/X+NE+IsX5WiVpEUiuDLSS7DHPG6BKUnN0nhLddp5USey6C3zZydH4OP8I9hxpRE5GnNevH3CggxAiG8ANAMZ2P7QLwAtSygKv90SkALVvkvU3kNTqttNruLg/3XaskoJTdlI0AGBLcZ1PoeSythJCzATwNQA7gFcALEXXCLzVQojpvjSWSC+eVElqBJKn1Bj67WmgaNVtxyopOCTF2BBlC8GeI40+vX6gSuk+AD+SUn7e67H/CSE+A/AAuu5bItKM2t12brfjQyApcZNsoHbbUWASQiA9LgIFFU0+vX6gq1DZfQIJACCl/BLASJ/2RuQjva8j+XIvUqB123kaSJ5223Exv8CVHh+B/ZXKh9JAtVezT3sj8oHRA8mfbruBqB1Irqokf4aRe9pt5w677swtPS4cFY3taPBhZoeBuu+yhBCP9/O4AJDh9Z6IDEitQPK3SvI3kNwFi7eB5Gu3nSuskgJbenwEAOBAZTMmZ8V79dqBQumeAb72O6/2QuQjNaskLQPJSJQY2OBptx2rpOCUFNP1PSypbVEulKSUi3t/LoSwSSk55x1pxoyB5IqRqiRvXqdGILFKCnyDoru+90fqvb8W6/Z2WyHEDCHEdgD7uj+fJIT4l9d7IvKCWQPJ6EtSeNNtp0cgsUoKDFFhIQi3WnC4ToVQAvAEutZSqgYAKeVWAKd4vSciDwVSILl8rg5VklaB5CsGUuAQQmBQtA1l9a1ev9aTULJIKQ/1eczh9Z6IPGDWQHLFSFVSf9RYRZbddgQAiVFhOFynTigVCyFmAJBCiBAhxK8B7PV6T0RuqDmFkNqBpESV5C9fqiRPsduOvBUbbkVti/dDwj0JpZsB3AZgCIByALO6HyNSjJpLUegVSEapkox+HYkCU5TNirrWDq9f58kifxVSyiuklIO7P66QUla5e50Q4kUhREX3qrV9v3a7EEJ2r2JLQS4QA0lr3o640yqQPMEqKTBF26xobLXD6ZRevc7tcuhCiOfRtX7SMaSUN7p56RIAT6JrMtfe28tC11pNRR63kgKSv911Rg4kV1WS1hOv9lclaRlI7LYLXlE2KySAxjY74iJDPX6d21AC8Fmv/w8HcCGAYncvklKuFkIM6+dLfwdwJ4B3Pdg3BSi1A8nt61Wa7VtrZg4kCmxRthAAQH1rp7KhJKV8rffnQoilAFZ62b6e154PoFRKubXP0ur9PfdGADcCQGbWEF92Rwal91LmSgSSEbrtlBjYoNbQb08CiVXSsXqf81IzMnVujf9CLF1Xh+xOp1ev836tWmA4gKHevkgIEQng9wDu9eT5UsrnpJTTpJTTBg3mpadAEQyB5M8AB0+vDfkya0Pf1/gbSP4MbGAgHa/3OS8h0fznvJ6yw7srSp5dU6rttV0LgBoAd3u5H6BruYvhAHqqpEwA3wshZkgpj/iwPTKZYAgkJdiycl0Gmydh5Gm3XV8MJFJST2eYlAoOdBBd6TEJQGn3Q07p7R5wtGHbAST32nYhgGmejOQjc9Pq+pGageRpGCk1DNzX+et8vY7EQCKlie5aydvEGLD7rjuA3pZSOro/PN68EGI5gPUAxgghSoQQi7xrGgWCYAokvTGQyFhkr/96zpPRdxuFEFOklN971Rwpr3Tz9WHebE8r7k6i/MXyHAPpeOGZ41UZFu7PSLu+GEikhJbOrtnoYsI9iZkfuHy2EMIqpbQDmAvgBiHEfnStOCvQVURN8bm1BuTpCbTnefwlG5gW148AcwWSGpSerYGBREppbu8Kpdhwz4eDAwNXShsBTAFwgc+tMglfTqDby9v5y+aC2QNJ7TBSqlpiIJGR1bZ0ICIsBJFhIV69bqBQEgAgpdzvT8OMzp8TKKum4zGQPONPMA10DxIDiYyisrEdWQkRcHdPal8DhVKSEOI2V1+UUj7u1Z4MSKlZqVk1dWEgead3uAwUUJ7eCMtAIiM50tCGcamxXr9uoFAKARCNH+6BCihqLpMQjIwQSAPRKpAGusdoIP4sLeFq+DgDifTS1ulAWV0rLpni/cwUA4VSmZTyQd+bZVxqBFIwV0tGCSR3E6y6YoYBDa74E0iuMJDIX4XVzXBKYGJmnNevHeg+JVZIBtq2UZkhkLSeYNXXG1+93Ye/FZLSayIxkKjHjsMNEADyhiR4/dqBKqXTfG5REAumiikQAkmtKql3YGg1y4Ore5DYZUda21ZSh4mZcUiMCvP6tS5DSUpZ41erDMibk2jfkyan4T+WVlWhWl12gHbddn3DxJOQ8rba4uSqZBS1LR0oqGjCL08d5dPrvbvVNsAN9Nd776+5+0UN9GrJCDM1eMKo6yIp2b3naXUEMJBIG+v3V8MpgfMmpfv0+qAJJXcnUm+6k/LL24I2mIwUSP5USWbmbqogBhLpRUqJ1fsqkZsRh+zkaJ+2ERShpGQg9X5NsHXpmSmQjHZPkr88mbPO1eg6LmFOWtlzpBGHqlvwxwvd/7y6EhShNBB/ltV2F0yBVC0F48hCPXkSQj0GGurNQCItfZx/BHERobgoz/eVc4M6lPwJpN7bCOSKSakwMlKVZATehI4r3oYRwO46Uk9JbQu+LazBz08ZiQgv57vrLeBDydVJ1ZOTZM8JcqC/LHu25eqX2szVUiBXR+FDclTrwlMicFxxdwOsqzACGEikrje/L0FEWAgWzR3h13YCPpS81d9f6r0fc/WLHUgVk9JhFOiDG9QMIcDzmRi8rY4ABhIpo7C6GRsO1OAXp2T7dG9SbwEdSt5WSZ6cFHdUtLutnPprh1l+sQO5OlKKnpVQXwNVRgADidQnpcTS9YcQHxmKG+b5VyUBAR5K/ekvkLy9H8ZVt56ZqyW1wkire5J84U0XnhpB5G0A9eZPGAEMJFLOpkO12FnWgIcWTkBcpHcL+vUn6EJJSb5UTUZklupIja673mHTE1BGC6Ae7oKoB6sj0kprhwOvrC/EqORoXDljiCLbDNhQ6u9Eq0SV5I6rasmoXXhqBpISoxu1pGQYaRlCPVgdkdZe31SM6qYOPHvNNFhDBprf23MBG0qe6C+Q3P013vdEYeZqySwVkln4G0TehhDgPogAz+dtZCCRNwoqGvHJjiO45oShmDrU+9nAXQnaUPIlkHqe4y6YzHxtSSlaLNznSlpammb3Kql5XcgVb/4IYiCRGuxOJ55fcxDJsTbcceYYRbcdkKHkadfdMV/34npFf8HkabuM8svPKsl3vgaRryHUQ40wAhhI5L33t5ahqKYFz187DTHh/g9u6C0gQ8mdvn+Z+3IBvW8wsVpSX05yuMffq57gULJi8iWM/AkiX7qFGUaktoNVzXjr+xKck5uGM8anKL79oAgls11wV5vaVZKR3m9/u/LMEESA9+t9MZDIFx12J576sgCJUWF4+AJ17tcLilAaSH9/ebs6ifU9QfnSjWekLjyjmJBsU/U+pd7fN3cBpXXXnD+DZHypxPmzR/54dWMRSutasXTRDCT4OXODK0EXSu5OfgOdtMrKygY8aZl5JJ5S1K6SvOnC648SQ7V7U2vEnMv9+dElzEAif2wprsMnO47gp3OGY96oJNX2E3Ch5NWS531Obp508/QNpoGqJSNeVwqEAQ7+BpMS+/eGv3+o+PszxDAif9W3duLZ1fsxKjkady5QdrRdXwEXSr7y5rqDu4qJ1KdHMGkZRkr8McMwIiU4nRJPfVmAlnYH/nl9HsJDfV+WwhMBH0paXHTvXS2xC883vlxX0iKYzBZEAMOIlPXW5hJsL63HoxfnYnx6rOr7C/hQ6q33Sa/3ycxVldRWlO9y6hl/qiUOdlBOT2goGU4MIqIu20rq8Nb3pbh4SgYum5alyT6DKpS80TM5p7+TdBrxupJR+TMKr3eQeBtQWo+cU/LngWFEaqluasdTXxZgdEoMHr4gF0IITfbLUOpHf8sZDFQ1kXKUGB7u78wJ7qh9U6s7DCJSW6fDiSe+2Ae7U+Lpq6f4tby5t5SZ1tUgfBlZ1rfrbqD1dfp+rfdrjbgiqh6UOPka8ZrchGTb0Q9P5aSEH/1QQm6KjYFEqpNS4vk1B7C3vAl/u3QSRiZFa7r/oKyU1A4QDnbwX8/7p8fif/21wxusisjMPthWhjX7qvCb00fjrFztRxkHZSiRunJSwhUb9ahHODGIKFh9d6gWyzcW4dyJabj1tGxd2sBQ6sWTpbF5bckzPSdps4ST3oMWGESkt6KaFjz55T7kZMThr5dM0mxgQ18MJVKVklUTcGx4+BNQes+y0INhREZQ29KBxz7dg5jwUDx/7TRNBzb0xVAKMrkpNs2nGlK6auqh9XU7BhEFopYOOx5dsRtN7Xa8fuMJSI3T9xYWhhJpRq1wUhun/KFA1elw4vGVe1Fa24rFP56O3Mw4vZvEUCLtmSWcGEYUyJxS4umvCrDjcAP+fvkknDRavZm/vcFQIt30PukbJaAYRBQMpJRYuv4QNhyowe/OHosL8zL1btJRDKVewofkuB2BFwgj7/S4ruSOngHFa0UUbN7fVoYVO45g0dzhuGHeCL2bcwyGEhmOFgHF+4ooWK3eW4nlG4tw3sQ0/P7scboN/XYlKENJ70XijMCI1VJ/+gsPb4NKrQlxGUZkNluK6/DcmgOYPXIQHrtsEiwWYwUSEGCh5MuJNi0tzasF/jzBKYbUpees6wwiMquCiib847O9GJMSg2evmQqbVb97kQYSUBOyKmGga0Z9v9Z7PSW1Z6ZWA0+wnuNkqGRmZfWt+Osnu5EUY8OSn05HTHio3k1yiaHUj/6CKRAGOPSHJ9qBMYzI7OpbO/Hoit0IsQgsXTQTyTHG/gM6oLrvlKRUCPXX1WS0k5xZri9pyWjfIyJfdNid+Nune1Df2onlN8zC8MFRejfJLdUqJSHEi0KICiFEfq/H/iqE2C2E2CaEeFsIEa/W/vvT+1pP7+42X5Y193UpdKPiSfiHqojvBQUCp5R46quCrmtJl+chb0iC3k3yiJrdd0sALOjz2EoAOVLKiQD2ArhHxf0D0OaiuBmvJ/UnGE/GDCIKVMs3FmHjwRr8/pxxWJCTqndzPKZa952UcrUQYlifxz7t9ekGAJeotX9veTMKb6Aqyewj73pOzv5253l6kter25AhRIHs051H8MG2Mlx7wlAsmjtc7+Z4Rc9rSj8F8JqrLwohbgRwIwBkZg3xeKPeXB/pe7+SJ8EUaN12rvQ9abt6T/09ubt6vdJhxRAio+t9zoTH5gYAAA09SURBVEvN8H3an52H6/HyukKcOjYZ95473nA3x7qjSygJIX4PwA5gmavnSCmfA/AcAEyeMlUqte8JyTZFF4obqOtOz/tplKb1Sd3bSouhQ2bX+5w3fmKeT+e8muYOPPFFAYYNisITV+bBGmK+Adaah5IQ4joA5wI4TUqpWNgopacS6lsxeVIhedJ1x5Onsvh+EnWxO5144ot96HA48ew1UxFtM+fgak1bLYRYAOAuACdJKVu02u9Aq5+6mnLIkxAKlAEORGR+r35ThD1HGvHElXkYlRKjd3N8puaQ8OUA1gMYI4QoEUIsAvAkgBgAK4UQW4QQz6i1/4H0rWgYLkRkZuv3V+Pj/CP48exhOH9Sut7N8Yuao++u7OfhxWrtT2t9g+y4oAug60lEZFyVjW14fs0BTBkSj9+dPU7v5vjNfFfBPNDfdYa+IeFPteRrZcXrH0SkJKeU+Peq/bBYgH9ekYcwq/lP6eY/AgV5Ejb9Pcfs9yYRkTmtyD+CXWWNuO+8CchKjNS7OYow5/AMhfQ3PLx36PQMgPC2MmLXHRGprbKxDa9vKsapY5Jw6VTjLGfur6CqlLwNi5zkcLeB5GmVxK47IlKKlBIvrS2ExSLw0IW5prtBdiABG0qehoA/XW/9vZZVEhGp7buiWmwursNvzxiNjPgIvZujqIANJVf6Cw1fgonXkYhID3aHE8s2FGFkUhSumz1M7+YoLqBDyZsuM29CxtVzXVVJ7LojIqV8urMcRxracO95ExBqwmmE3Am8I/KAq/DwJJhYIRGRXto6HXh3aynmZA/CSaOT9G6OKoJ69F1/fA0dXksiIrWt2HEEDa123D5/jN5NUU3AV0quus6UDJGBtsWuOyJSQrvdgY+3l+GUMUmmWUXWFwEfSgNRIphYIRGRFlbtrURDmx03n5ytd1NUFRShNFC14k+ouHstqyQiUoJTSny8/QgmZ8Vj+rDArZKAIAkld3wJJgYSEWklv7QeRxra8JM5wwLqRtn+BE0ouQuJnJRwj8LJ0+cRESnl890ViI8MxYKcVL2bojqOvuujd+D0LAzobQixSiIipTS12/H9oVpce8Iw2KwhejdHdUFTKQHeh4UvVREDiYiUtPFgDexOiQvyzL14n6eCKpQAhgYRmcuGA9UYOigSuRlxejdFE0EXSoB6wcTAIyIltXTYsbOsAQtyUgN+gEOPoAwlQPkAYSARkdK2ldTD4ZQ4fVyK3k3RTNCGEqBckDCQiEgNOw7XI9pmRV5WvN5N0UxQhxLgX6DkptgYSESkmr3ljZg6NAHWAJwN3JXgOdIB+BIsDCMiUlNTux1FNa2YNjSwZ3Doi/cpdesJme3l7R49j4hITXvLGwEA04Yl6twSbTGU+mDoEJER7DnSCKtFYHIQXU8C2H1HRGRIRTUtGJUSjYiwwJ/FoTeGEhGRAR2ua0V2cozezdAcQ4mIyGCcUqKysR3ZSdF6N0VzDCUiIoPpdDghAWQnM5SIiEhnHXYnAIYSEREZgN0pAQCZCRE6t0R7DCUiIoNxOCWibVZE2YLvrh2GEhGRwdgdEimxwXnPJEOJiMhgHE6J1DjvFhgNFAwlIiKDsTudSIlhKBERkQE4pERiVJjezdAFQ4mIyGCkBGIjQvVuhi4YSkREBhQbHnwj7wCGEhGRIcVFslIiIiKDiLExlIiIyCAig2zJih4MJSIiA7KFBufpOTiPmojI4GxWVkpERGQQ4ayUiIjIKFgpERGRYfCaEhERGQYrJSIiMozQEKF3E3TBUCIiMiCLYCgREZFBMJSIiMgwLMGZSeqFkhDiRSFEhRAiv9djiUKIlUKIfd3/Jqi1fyIiM2OlpLwlABb0eexuAJ9LKUcB+Lz7cyIi6iNIM0m9UJJSrgZQ0+fhhQBe7v7/lwFcoNb+iYjMTARpKml9TSlFSlkGAN3/Jrt6ohDiRiHEJiHEpuqqKs0aSESkh97nPL3boifDDnSQUj4npZwmpZw2aPBgvZtDRKSq3ue84KyRumgdSuVCiDQA6P63QuP9ExGRgWkdSu8BuK77/68D8K7G+yciMr4gLpXUHBK+HMB6AGOEECVCiEUAHgFwhhBiH4Azuj8nIqJeRBCnklWtDUspr3TxpdPU2icREZmbYQc6EBEFq+CtkxhKRERkIAwlIiKjCeJSiaFERGQwQZxJDCUiIjIOhhIRERkGQ4mIyGCC+T4lhhIRkdEEbyYxlIiIDEfq3QD9MJSIiAxGBnEqMZSIiMgwGEpERGQYDCUiIjIMhhIRERkGQ4mIiAyDoURERIbBUCIiIsNgKBERGUzw3qXEUCIiIgNhKBERGU0Ql0oMJSIig3HK4E0lhhIRkcFIADJIg4mhRERkQB0Op95N0AVDiYjIgNo6gzOUhBlKRCFEJYBDerfDA4MBVOndCJ0E67HzuIOLv8ddJaVc4O5JQogVnjwvEJkilMxCCLFJSjlN73boIViPnccdXIL1uLXE7jsiIjIMhhIRERkGQ0lZz+ndAB0F67HzuINLsB63ZnhNiYiIDIOVEhERGQZDiYiIDIOhpBAhxG+EEDuEEPlCiOVCiHC926QGIcSLQogKIUR+r8cShRArhRD7uv9N0LONanBx3H8VQuwWQmwTQrwthIjXs41q6e/Ye33tdiGEFEIM1qNtanJ13EKIXwoh9nT/vv9Fr/YFKoaSAoQQGQBuBTBNSpkDIATAFfq2SjVLAPS9qe9uAJ9LKUcB+Lz780CzBMcf90oAOVLKiQD2ArhH60ZpZAmOP3YIIbIAnAGgSOsGaWQJ+hy3EOIUAAsBTJRSTgDwmA7tCmgMJeVYAUQIIawAIgEc1rk9qpBSrgZQ0+fhhQBe7v7/lwFcoGmjNNDfcUspP5VS2rs/3QAgU/OGacDF9xwA/g7gTgToQgsujvtmAI9IKdu7n1OhecMCHENJAVLKUnT9xVQEoAxAvZTyU31bpakUKWUZAHT/m6xze/TwUwAf690IrQghzgdQKqXcqndbNDYawDwhxDdCiFVCiOl6NyjQMJQU0H0NZSGA4QDSAUQJIa7Wt1WkFSHE7wHYASzTuy1aEEJEAvg9gHv1bosOrAASAMwCcAeA14UQQt8mBRaGkjJOB3BQSlkppewE8BaA2Tq3SUvlQog0AOj+N2i6NIQQ1wE4F8CPZPDc9DcSXX+AbRVCFKKr2/J7IUSqrq3SRgmAt2SXjQCc6JqklRTCUFJGEYBZQojI7r+aTgOwS+c2aek9ANd1//91AN7VsS2aEUIsAHAXgPOllC16t0crUsrtUspkKeUwKeUwdJ2op0gpj+jcNC28A+BUABBCjAYQhuCcLV01DCUFSCm/AfA/AN8D2I6u9zUgpyMRQiwHsB7AGCFEiRBiEYBHAJwhhNiHrtFYj+jZRjW4OO4nAcQAWCmE2CKEeEbXRqrExbEHPBfH/SKAEd3DxP8L4LogqpA1wWmGiIjIMFgpERGRYTCUiIjIMBhKRERkGAwlIiIyDIYSEREZBkOJDEcI0eTBc+Z1z9K8RQgRoXJ7hvXMFC2EOFkIUS+E2Nw9U/RqIcS5vZ77MyHEtQNs62QhRDDdWE3kFaveDSDy0Y8APCalfMmTJ3ff1CyklM5ej4VIKR0+7HuNlPLc7m1MBvCOEKJVSvm5lNLdvUonA2gCsM6H/RIFPFZKZFjdVcVXQoj/da9btEx0uR7AZQDuFUIs637uHUKIb7vXNnqg+7FhQohdQoin0XVjc5YQokkI8aAQ4hsAJwghpnZPrPmdEOKTXtMlTRVCbBVCrAdwi6s2Sim3AHgQwC+6X3e/EOL27v+/VQixs7tN/xVCDAPwMwC/6a7w5qnzzhGZFyslMro8ABPQtRTIWgBzpJQvCCHmAvhASvk/IcR8AKMAzAAgALwnhDgR/9/eHYNWEURRGP6PmsqAr7EUtAohIlim0FIQwTSCjU1SClrbxE4sLBTBwkYLtVMQrBREHgqKhRExhaUiCNpEjGKEcCxmiueyJsZqjOeDB7vLnWUXHlxmZrm3lH+aAGZtnwSQtB14bfuspDFgCMzY/iTpOHCOUvH7OnDK9lDShXWe8QWlOGfXGWCP7RVJA9tLterDsu304YnokaQUrXtu+z2ApJfAbuBJJ+ZQ/S3U83FKknoHvLX9bCR2FbhTjyeAvZQyQVCaM36QtAMY2B7WuBvA4TWe8XdVol8BtyTdpdRMi4h1JClF61ZGjlfp/88KOG/76i8Xy3LZ107s95F9JAGLtqc74wZsrHHdfvoL8B4BDgJHgXlJUxu4Z8R/KXtKsRncB+YkjUNpTy/pTxoNvgF2Spqu48YkTdleAj7XJUIoH1X0krQPmAeudK5vAXbZfkTpzjqgzOC+UIq4RkSPzJTin2f7gaRJ4GldhlsGTlBmVmuN+yHpGHC5LtltAy4Bi8AscE3SN0rSG3VA0gKl7f1H4LTth52YrcDNel8BF+ue0j3gtqQZyp7V479/84jNJ1XCIyKiGVm+i4iIZiQpRUREM5KUIiKiGUlKERHRjCSliIhoRpJSREQ0I0kpIiKa8RPPUaElcWhx2wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.jointplot(np.log(tf[\"InferredDist\"]), np.log(tf[\"TrueDist\"]), kind=\"kde\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now iterate over several threshold to see if one could differentiate, qualitatively, adjacent contigs against non-adjacent contigs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1029,\n",
       " 1162,\n",
       " 1020,\n",
       " 61,\n",
       " 0.4674610449129239,\n",
       " 0.9440366972477064,\n",
       " 0.6696210268948656,\n",
       " [Pair            idcChr1.ctg43;idcChr1.ctg48\n",
       "  InferredDist                         623487\n",
       "  TrueDist                             956309\n",
       "  Name: 6, dtype: object, Pair            idcChr1.ctg47;idcChr1.ctg49\n",
       "  InferredDist                         571740\n",
       "  TrueDist                             630265\n",
       "  Name: 135, dtype: object, Pair            jpcChr1.ctg451;jpcChr1.ctg460\n",
       "  InferredDist                           623487\n",
       "  TrueDist                               982109\n",
       "  Name: 251, dtype: object, Pair            idcChr1.ctg330;idcChr1.ctg338\n",
       "  InferredDist                           679917\n",
       "  TrueDist                               828578\n",
       "  Name: 328, dtype: object, Pair            idcChr1.ctg110;idcChr1.ctg119\n",
       "  InferredDist                           547500\n",
       "  TrueDist                               983402\n",
       "  Name: 473, dtype: object, Pair            idcChr1.ctg37;idcChr1.ctg47\n",
       "  InferredDist                        1004119\n",
       "  TrueDist                             988917\n",
       "  Name: 530, dtype: object, Pair            idcChr1.ctg338;idcChr1.ctg349\n",
       "  InferredDist                           679917\n",
       "  TrueDist                               783337\n",
       "  Name: 552, dtype: object, Pair            idcChr1.ctg411;idcChr1.ctg413\n",
       "  InferredDist                           547500\n",
       "  TrueDist                               598131\n",
       "  Name: 680, dtype: object, Pair            idcChr1.ctg40;idcChr1.ctg47\n",
       "  InferredDist                        1004119\n",
       "  TrueDist                             851616\n",
       "  Name: 692, dtype: object, Pair            idcChr1.ctg71;idcChr1.ctg80\n",
       "  InferredDist                         623487\n",
       "  TrueDist                             904514\n",
       "  Name: 859, dtype: object])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def CrossValidation(df, threshold=0, truethreshold=1000000):\n",
    "    TP, TN, FP, FN = 0, 0, 0, 0\n",
    "    TPS, TNS, FPS, FNS = [], [], [], []\n",
    "    for i, row in df.iterrows():\n",
    "        inferred_dist, true_dist = row[\"InferredDist\"], row[\"TrueDist\"]\n",
    "        if inferred_dist < threshold: # positive\n",
    "            if true_dist < truethreshold:\n",
    "                TP += 1\n",
    "            else:\n",
    "                FP += 1\n",
    "        else:  # negative\n",
    "            if true_dist >= truethreshold:\n",
    "                TN += 1\n",
    "            else:\n",
    "                FN += 1\n",
    "                FNS.append(row)\n",
    "    FPR = FP * 1. / (FP + TN)\n",
    "    TPR = TP * 1. / (TP + FN)\n",
    "    ACC = (TP + TN) * 1. / (TP + TN + FP + FN)\n",
    "    return TP, TN, FP, FN, FPR, TPR, ACC, FNS[:10]\n",
    "\n",
    "CrossValidation(tf, threshold=500000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 (0, 2182, 0, 1090, 0.0, 0.0, 0.6668704156479217)\n",
      "500000 (1029, 1162, 1020, 61, 0.4674610449129239, 0.9440366972477064, 0.6696210268948656)\n",
      "1000000 (1082, 594, 1588, 8, 0.7277726856095326, 0.9926605504587156, 0.5122249388753056)\n",
      "1500000 (1088, 292, 1890, 2, 0.8661778185151238, 0.998165137614679, 0.421760391198044)\n",
      "2000000 (1088, 160, 2022, 2, 0.9266727772685609, 0.998165137614679, 0.38141809290953543)\n",
      "2500000 (1089, 73, 2109, 1, 0.9665444546287809, 0.9990825688073395, 0.35513447432762835)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Threshold</th>\n",
       "      <th>TP</th>\n",
       "      <th>TN</th>\n",
       "      <th>FP</th>\n",
       "      <th>FN</th>\n",
       "      <th>FPR</th>\n",
       "      <th>TPR</th>\n",
       "      <th>ACC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2182</td>\n",
       "      <td>0</td>\n",
       "      <td>1090</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.666870</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>100000</td>\n",
       "      <td>691</td>\n",
       "      <td>1982</td>\n",
       "      <td>200</td>\n",
       "      <td>399</td>\n",
       "      <td>0.091659</td>\n",
       "      <td>0.633945</td>\n",
       "      <td>0.816932</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>200000</td>\n",
       "      <td>859</td>\n",
       "      <td>1764</td>\n",
       "      <td>418</td>\n",
       "      <td>231</td>\n",
       "      <td>0.191567</td>\n",
       "      <td>0.788073</td>\n",
       "      <td>0.801650</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>300000</td>\n",
       "      <td>952</td>\n",
       "      <td>1556</td>\n",
       "      <td>626</td>\n",
       "      <td>138</td>\n",
       "      <td>0.286893</td>\n",
       "      <td>0.873394</td>\n",
       "      <td>0.766504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>400000</td>\n",
       "      <td>997</td>\n",
       "      <td>1339</td>\n",
       "      <td>843</td>\n",
       "      <td>93</td>\n",
       "      <td>0.386343</td>\n",
       "      <td>0.914679</td>\n",
       "      <td>0.713936</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Threshold   TP    TN   FP    FN       FPR       TPR       ACC\n",
       "0          0    0  2182    0  1090  0.000000  0.000000  0.666870\n",
       "1     100000  691  1982  200   399  0.091659  0.633945  0.816932\n",
       "2     200000  859  1764  418   231  0.191567  0.788073  0.801650\n",
       "3     300000  952  1556  626   138  0.286893  0.873394  0.766504\n",
       "4     400000  997  1339  843    93  0.386343  0.914679  0.713936"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roc = []\n",
    "for i in range(0, 3000000, 100000):\n",
    "    c = CrossValidation(tf, threshold=i)\n",
    "    if i % 500000 == 0: \n",
    "        print i, c[:-1]\n",
    "    roc.append(tuple([i] + list(c[:-1])))\n",
    "    \n",
    "rf = pd.DataFrame(roc, columns=[\"Threshold\", \"TP\", \"TN\", \"FP\", \"FN\", \"FPR\", \"TPR\", \"ACC\"])\n",
    "rf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'ROC curve')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8XXWd//HXO0nTdEk3WtrSFboAZW9DRXBBUWQZQR0UUFTGhcFxGcVZcFzGQcdxHZURRTZZfoOAjqOVQVERF5ACSYFCgUKBpk1p6Xabpkv2z++PcxrSkKS3y83Nvff9fDzuI2f53nM+p0nv536/33O+X0UEZmZmAGX5DsDMzAYPJwUzM+vipGBmZl2cFMzMrIuTgpmZdXFSMDOzLk4KZmbWxUnBio6klZJ2StomaZ2kGyWN7FHmZEm/l9QkqVHSLyXN61FmlKTvSFqVHmtFuj5+YK/IbOA4KVixemtEjASOB04APrNrh6RXA78BfgEcAhwKPAbcL+mwtEwlcA9wFHAGMAo4GdgELMxV0JIqcnVss2w4KVhRi4h1wN0kyWGXrwM3R8R3I6IpIjZHxOeAxcAX0zLvA6YDb4+IJyOiMyLWR8SXIuKu3s4l6ShJv5W0WdJLkv4l3X6jpC93K3eqpIZu6ysl/bOkpcB2SZ+T9NMex/6upCvT5dGSrpe0VtIaSV+WVL6f/1RmgJOCFTlJU4EzgRXp+nCSb/w/6aX4HcCb0+U3Ab+OiG1Znqca+B3wa5Lax2ySmka2LgTOBsYAtwBnSRqVHrsceBdwa1r2JqA9PccJwOnAh/biXGZ9clKwYvVzSU3AamA98K/p9nEkf/dre3nPWmBXf8FBfZTpy18B6yLiWxHRnNZAHtyL918ZEasjYmdE1ANLgLel+94I7IiIxZImkiS5T0bE9ohYD3wbuGAvzmXWJycFK1Zvi4hq4FTgCF7+sM8AncDkXt4zGdiYLm/qo0xfpgHP7VOkidU91m8lqT0AvJuXawkzgCHAWklbJG0BfggcvB/nNuvipGBFLSL+CNwIfDNd3w48ALyzl+Lv4uUmn98Bb5E0IstTrQZm9bFvOzC82/qk3kLtsf4T4NS0+evtvJwUVgMtwPiIGJO+RkXEUVnGadYvJwUrBd8B3ixpV2fz5cD7JX1CUrWksWlH8KuBf0vL3ELyAfw/ko6QVCbpIEn/IumsXs5xJzBJ0iclDU2P+6p036MkfQTjJE0CPrmngCNiA/AH4EfACxHxVLp9LcmdU99Kb5ktkzRL0uv34d/F7BWcFKzopR+wNwOfT9fvA94CvIOk36CepMP2NRHxbFqmhaSz+Wngt8BW4CGSZqhX9BVERBNJJ/VbgXXAs8Ab0t23kNzyupLkA/32LEO/NY3h1h7b3wdUAk+SNIf9lL1r6jLrkzzJjpmZ7eKagpmZdXFSMDOzLk4KZmbWxUnBzMy6FNzgW+PHj4+ZM2fmOwwzs4JSV1e3MSIm7KlcwSWFmTNnUltbm+8wzMwKiqT6bMq5+cjMzLo4KZiZWRcnBTMz6+KkYGZmXZwUzMysS86SgqQbJK2X9EQf+yXpynQy9KWS5ucqFjMzy04uawo3kkx43pczgTnp6xLgBzmMxczMspCz5xQi4k+SZvZT5FySydMDWCxpjKTJ6XjxZmZZ6egMmts62NnWQXNbB20dhTHyc0Swo7WDrc1tNDW309m5e9wBbGtpZ0NTCy1tHQCcduREjps2Jqdx5fPhtSnsPgVhQ7rtFUlB0iUktQmmT58+IMGZ2YHX0RnsbOtge0s7Tc3tbGtp71re3pKsd73SbU1pme0t7ekHf2dXAiikJLA/pOTnwaOqijopqJdtvf52I+Ia4BqAmpqa4v8LMMuzjs5gR2s721s62N7azo5dP9NtXfta2tne2rH79tYOmls72NHWzo7WDna2Jt/id7R20NremdX5K8rEyKoKRg59+TV6eCWTh5RTNaSMYZXlDK0oZ1hlOVUV5QyrLKNqSLI8pEKo14+XwWdYZTmjqoZQXVVBRfkrYx5RWcGE6qFUDSkfsJjymRQaSCY732Uq8GKeYjEreK3tnbywcTvLX2pi/dZmWto7u75NN7ely922tbR10tz+8s/u5Vqy/PAGqKwoY0RlOcMrKxgxtJxhlRUMH1LOwdVVDKssZ/iQ5MN7WGU5w4aUM7yynBFDd//AHzG0guqql5eHVpQhFcYHe7HJZ1JYBHxM0m3Aq4BG9yeY7VlHZ7Bq8w6Wr2vimZeaWP5SE8+sa+KFjdtp79EuLUFVRfLtumpIOVVDyhlaUdb1c+zwypf3dSs3rLKckUMruj7oh1dWMCL9MH95vYJhleVUVvjO9mKSs6Qg6cfAqcB4SQ3AvwJDACLiauAu4CxgBbAD+JtcxWJWiCKCFxubeWZd+sGfvp59adtu3+SnjxvO3InVnH7UROZOrGbuxGoOGTOMqiFlVJb7G7ftnVzefXThHvYH8NFcnd+skEQE9Zt2sHRNI0tXb2HpmkaeenErTS3tXWUmjapi7qRq3vfqg5gzsZrDJ1Yz++CRjBhacIMd2yDmvyazARYRrG1sZmnDFpY2NKavLWxtThLA0Ioy5h0yirfPn8LcidUcPqmauQdXM3r4kDxHbqXAScEsxzZua2FpwxYeW93I42uSJLBxWwuQ3GVzxORqzj72EI6bOppjpo5m7sRqhpS7nd7yw0nB7ABq3NHG42saeaxhC4+nNYAXG5uBpNN3zsEjef3cCRw3bTTHTBnNkZNHDejthmZ74qRgto+2t7TzxJrGNAk08njDFlZu2tG1f+ZBw6mZOY5jp47m2KljOOqQUW7/t0HPf6FmWejoDJ5au5UlqzI8tjqpAazYsI1I7wA9ZHQVx04dw7tOnMaxU8ZwzJTR7gOwguSkYNaLbS3tPLIqQ+3KDHX1GR5ZlWF7azL+zPiRlRw7dQxnHzuZY6eO5pgpY5hQPTTPEZsdGE4KVvIigjVbdlJXnySA2pUZnl63lc6AMsHhk0bxjvlTqZk5lvnTxzJ17DDf+29Fy0nBSk57RydPrW2itn4ztfUZ6lZmWLc16QweUVnOCdPH8vE3zqFm5liOnzaG6io3A1npcFKwore1uY0l9RmW1Georc/w6Oot7EibgqaMGcaJh46jZsZYFswYyxGTqqnw7aBWwpwUrKhEBA2ZnUktIO0PWP5SExFQXiaOnFzNu2qmsWDGWGpmjmXy6GH5DtlsUHFSsIIWEbywcTv3r9jIA89vonZlhvVNyYNh1UMrOGHGWM46ZjILZiRNQb4l1Kx//h9iBWf91mbuf24j96/YxP0rNrI2fThsyphhnDzrIBbMTJqD5k6sprzMHcJme8NJwQa9puY2Hnx+M/et2MhfntvIMy9tA2DM8CGcMms8J88+iNfMHs/0ccN9V5DZfnJSsEGnpb2DJfVb+MtzG7l/xUYea2ikozOoGlLGwkMP4q/nT+WU2eOZN3kUZa4JmB1QTgqWd52dwZNrt3L/io3ct2IjD6/cTHNbJ+Vl4tipo/m7U2dx8qzxzJ8xhqEVHifILJecFGzA7Zo74P60JvDAc5vI7GgDkgHjLjhxOqfMHs+rDhvHKD8jYDagnBRsQGxoaulqDrp/xSbWbNkJwOTRVZx25EROmX0QJ88az8RRVXmO1Ky0OSlYTmxraeehFzZ13SH09LomAEZVVXDyrPFc+vrDOGX2eA4dP8Kdw2aDiJOCHTCbt7dyywP1/PnZDTy6egvtnUFlRRknzhzLP51xOKfMGs/RU0b7NlGzQcxJwfZba3snNz+wkivveZamlnaOnTKaS16X1AQWzBjrSWTMCoiTgu2ziOCep9bz73c9xQsbt/O6uRP4/NlHMmdidb5DM7N95KRg++SptVv58v89yf0rNjFrwgh+9Dcn8obDD853WGa2n5wUbK9s3NbCt37zDLc/vIpRw4bwb+ccxbtfNd0TzZsVCScFy0pLewc33r+S7/1+BTvbOnj/yTP5+9PmMGZ4Zb5DM7MDyEnB+hUR3L1sHV+562lWbd7BaUcczL+cfSSzJozMd2hmlgNOCtanJ9Y08qU7n+TBFzZz+MRqbvngQl47Z0K+wzKzHHJSsFdYv7WZb/5mOT+pa2Ds8Eq+/LajueDEaZ6RzKwEOClYl+a2Dq6/7wW+f+8KWjs6+fBrD+Ojb5jN6GEef8isVDgpGBHBnUvX8tVfPc2aLTt5y1ET+cyZRzJz/Ih8h2ZmA8xJocQ9tnoLX7rzSWrrMxw5eRTfeOexnDxrfL7DMrM8cVIoUesam/n6r5/mZ4+sYfzIoXztr4/hvAXTPC6RWYnLaVKQdAbwXaAcuC4ivtpj/3TgJmBMWubyiLgrlzGVup2tHVzzp+e5+o/P0RHBR06dxd+dOotqz1tgZuQwKUgqB64C3gw0AA9LWhQRT3Yr9jngjoj4gaR5wF3AzFzFVMo6O4NFj73I1379NGsbmzn7mMlcfuYRTBs3PN+hmdkgksuawkJgRUQ8DyDpNuBcoHtSCGBUujwaeDGH8ZSsuvoMX7rzSR5dvYVjpozmuxecwMJDx+U7LDMbhHKZFKYAq7utNwCv6lHmi8BvJH0cGAG8qbcDSboEuARg+vTpBzzQYtWQ2cHXfr2cXz72IgdXD+Wb7zyOd5wwxZPdm1mfcpkUevvkiR7rFwI3RsS3JL0auEXS0RHRudubIq4BrgGoqanpeQzrISL4/h+e48p7ngXgE2+czd++fhYjhvq+AjPrX1afEpJGAZOBncDqiMjmg7kBmNZtfSqvbB76IHAGQEQ8IKkKGA+szyYue6WI4Io7n+RH96/krGMm8dmz5zFlzLB8h2VmBaLPpCCpGvgI8G5gJLARqAIOknQf8P2I+HM/x34YmCPpUGANcEF6rO5WAacBN0o6Mj3+hn28lpLX2Rl88ZfLuPmBej5wyqF8/q+O9PzHZrZX+qsp/C/w38BpEbFp10YlnzILgfdKmhMRN/T25ohol/Qx4G6S201viIhlkq4AaiNiEfBp4FpJnyJpWro4y1qI9dDZGXz250/w44dW8bevO4zLzzzCCcHM9poK7TO4pqYmamtr8x3GoNLRGXzmZ0u5o7aBj75hFv9w+uFOCGa2G0l1EVGzp3J7PeylpNmSfrBvYdmB1tEZ/ONPHuOO2gb+/rQ5Tghmtl/6TAqSjpZ0l6RHJX1R0gRJtwN/Ap4fuBCtL+0dnXzq9kf52SNr+PSb5/KpN891QjCz/dJfn8J16esBkjuElgA/AWZFxM4BiM360dbRySdve5T/e3wt/3zGEXzk1Fn5DsnMikB/SaEqIq5Ll5dJ+iTwTxHRPgBxWT9a2zv5+I+XcPeyl/jc2Ufyodcelu+QzKxI9JsUJB3Dyw+hbQOOTO8+IiKW5jo4e6WW9g4++t9L+N1T6/niW+dx8SmH5jskMysi/SWFjcD3+1gP4HW5Csp619zWwaX/r44/LN/Al992NBedNCPfIZlZkekzKUTEawYyEOvfztYOLrmllvtWbOSr7ziGCxZ6DCgzO/D6u/tolqSfpncf3SJp8kAGZi/b0drOB258mPtWbOQb5x3nhGBmOdPfcwo/Au4B3kMy3PV/DUhEtpttLe1cfMPDPPjCJr79ruM5b8HUfIdkZkWsvz6FURGx6yG1ZZKWDERA9rKm5jYu/tHDPLp6C9+94ATeetwh+Q7JzIrc3tx9NKz7uu8+yq3GnW2874aHWLamke9deAJnHuPWOzPLvf6SwgZ891FebNnRynuvf4in123l+++Zz+lHTcp3SGZWIvpLCpdFxMMDFokBsHl7Kxdd9yAr1m/jh+9dwBuPmJjvkMyshPTX0fzDAYvCANi4rYV3X7uY5zZs49r31zghmNmA66+m4JHVBtD6pmbec+2DrM7s4IaLT+SU2ePzHZKZlaD+ksKhkn7W186IeEcO4ilJL21t5sJrF7OusZkfXbyQV886KN8hmVmJ2lNH81UDFUipenHLTt597WI2NLVw0wcWcuLMcfkOycxKWH9JoSki7hmwSEpQQ2YHF167mC3b27j5g69iwYyx+Q7JzEpcf0lh9YBFUYJWbUoSQlNzG7d86FUcP21MvkMyM+v37qP/6O+NkkZKmneA4ykJKzdu5/xrHmB7azu3fvgkJwQzGzT6qym8R9I3gF8BdSR9DFXAbOAN6c9/yHmERea5Ddt497WLaesIbv3QScw7ZFS+QzIz69Lf0NkflzQeeCfwXmAysBN4CrgpIv4wIBEWkWdfauLd1z1IRPDjD5/E4ZOq8x2Smdlu+qspEBEbgR+kL9sPT6/bynuufZCyMnHbJScx+2AnBDMbfPrrU7ADZNmLjVx4zWIqyp0QzGxw67emYPvv8YZGLrr+QUZUlnPrh09i5vgR+Q7JzKxPTgo59OjqLbz3+gcZVTWE2y45iWnjhuc7JDOzfu2x+UjSMEmfkXR1uj5b0pm5D62wdXYGH765lrHDK7n9b50QzKwwZNOncAPJ4HivSddfBL6Ss4iKxLPrt7GhqYVPnDaHqWOdEMysMGSTFOZExFeANoCI2IFHUN2juvoMADUeusLMCkg2SaFVUhXJbGtIOhRozWlURaCuPsNBIyqZcZBrCWZWOLJJCl8Cfg1MlXQTcC/wL9kcXNIZkpZLWiHp8j7KvEvSk5KWSbo168gHuSWrMsyfMRbJlSozKxx7vPsoIn4lqRY4maTZ6B8jYv2e3iepnGTo7TcDDcDDkhZFxJPdyswBPgOcEhEZSQfv43UMKhu3tfDCxu2cf+K0fIdiZrZXsrn76DcRsSEifhERP4+I9ZJ+k8WxFwIrIuL5iGgFbgPO7VHmw8BVEZEByCbZFIIl7k8wswLVZ1KQVClpFDBRUrWkUelrKjA9i2NPYffhtxvSbd3NBeZKul/SYkln9BHLJZJqJdVu2LAhi1PnV92qDEPKxdFTRuc7FDOzvdJf89FHgcuAg4FlvHzH0Vbg6iyO3VtjevRy/jnAqcBU4M+Sjo6ILbu9KeIa4BqAmpqanscYdJbUZzh6ymiqhpTnOxQzs73SZ00hIr4dEdOAf46I6RExLX0dFRHfyeLYDUD3RvWpJM849Czzi4hoi4gXgOUkSaJgtbR38FhDIwumu+nIzApPNh3N35F0BDCPZD6FXdv3dKfQw8Cc9BbWNcAFwLt7lPk5cCFwYzpM91zg+ezDH3yWvbiV1vZOamY6KZhZ4dljUpD0OeB04AjgbuAtwH1Av0khItolfSx9TzlwQ0Qsk3QFUBsRi9J9p0t6EuggubNp0/5cUL7t6mSe75qCmRWgbAbEOx84HlgSEe+VNBn4YTYHj4i7gLt6bPtCt+Ug6be4LOuIB7m6+gzTxg3j4FFVey5sZjbIZPPw2s6I6ADaJVUD64DDchtWYYoIausz7k8ws4KVTU3hEUljSAbGqyW5+2hJTqMqUA2ZnWxoamHBzHH5DsXMbJ/0mxSUjNHwxfQW0ask3Q2MiggnhV7sGgTPNQUzK1T9Nh+lbf53dltf4YTQt7r6DCMqyzl8kqfbNLPClE2fwkOS5uc8kiJQW5/hhOljKS/zIHhmVpiySQqvIUkMyyUtkfSIJNcWemhqbmP5uq0s8HhHZlbAsuloflvOoygCj61upDNwUjCzgpbNE83PDUQgha6uPoMEx08fk+9QzMz2WTbNR5aF2vrNHD6xmlFVQ/IdipnZPnNSOAA6OoNHV21x05GZFbyskoKkqZLekC4PlTQit2EVlmfXN9HU0u6kYGYFL5uZ1z4ALAKuSzfNAH6Ry6AKTddDa04KZlbgsqkpfAI4iWR4CyLiGZKJdyxVtzLD+JGVTB83PN+hmJntl2ySQnM6xzIAksrpfVa1klW3KsOCGWNJRgUxMytc2SSF+yX9E1CV9ivcTrehL0rdhqYW6jftcNORmRWFbJLCPwFNwNPA3wP3AJ/NZVCFZMkq9yeYWfHI5onms4DrIuIHuQ6mENXVZ6gsL+OoQ0bnOxQzs/2WTU3hXcAKST+S9Ja0T8FSdfUZjpk6mqoh/mcxs8K3x6QQEe8F5gK/BD4APC/p6lwHVgha2jt4vKHRTUdmVjSyaT4iIlok/QLYCZST1B4uzWVgheCJNVtp7ehkvifVMbMikc3Da2+SdB3wHHARcDMwKdeBFYK6+s2AO5nNrHhkU1O4FLgN+HhE7MxxPAWlrj7DjIOGM6F6aL5DMTM7ILIZOvu8gQik0EQEdfVbeN2c8fkOxczsgOkzKUj6Y0S8XlIGiO67SKZvHpfz6Aax1Zt3snFbC/PddGRmRaS/msIb0p/+KtyLWvcnmFkR6rOjOSI608XrI6Kj+wu4fmDCG7zq6jNUD61g7sTqfIdiZnbAZPPw2rHdV9KH107MTTiFo64+w/HTx1Be5kHwzKx49JkUJP1z2p9wrKTN6SsDbADuGrAIB6Gm5jaWv9TkpiMzKzr91RS+DkwAvp3+nACMj4hxEfGPAxHcYPXIqi1EuD/BzIpPfx3NsyPiWUm3AEft2rhrzoCIWJrj2AatuvoMZYLjp43JdyhmZgdUf0nhcuCDwFW97AvgdTmJqAAsWZXh8EmjqK4aku9QzMwOqP7uPvpg+vO1vbyySgiSzpC0XNIKSZf3U+48SSGpZu8vYWB1dAaPrNrCghmuJZhZ8clm7KN3SKpOly+XdIek47J4XzlJLeNMYB5woaR5vZSrJpkH+sG9DT4flq9rYltLu/sTzKwoZXNL6hcjoknSycBbSabj/GEW71sIrIiI59M5nm8Dzu2l3JdIOrWbs4w5r+p2zbQ2vaQf6DazIpVNUuhIf/4V8P2I+B8gmxHgpgCru603pNu6SDoBmBYR/c75LOkSSbWSajds2JDFqXNnSX2GCdVDmTZuWF7jMDPLhWySwlpJVwEXAHdJqszyfb091dU1hpKkMpLbXT+9pwNFxDURURMRNRMmTMji1LlTV59hwfSxXXdhmZkVk2yn4/wjcFZEZEjGQuqz07ibBmBat/WpwIvd1quBo4E/SFoJnAQsGsydzeubmlm1eYf7E8ysaGUzHec24EngVEmXAmMj4ldZHPthYI6kQ9PaxQXAom7HbYyI8RExMyJmAouBcyKidl8uZCAsqU/6EzwyqpkVq2zuPvoYcAcwPX3dIenv9vS+iGgHPgbcDTwF3BERyyRdIemc/Qs7P+rqM1RWlHH0lFH5DsXMLCeymXntEmBhWmNA0leAvwDf39MbI+IueoyTFBFf6KPsqVnEkld19RmOnTKaoRXl+Q7FzCwnsu0wbuu23kbvnchFrbmtgyfWbHV/gpkVtWxqCrcAiyX9D0kyeBtwU06jGoSeWNNIa0en+xPMrKhlM0fz1yXdC7w23XRpRDyc27AGn7q0k9k1BTMrZtnUFABa0ldn+rPk1NVnmHnQcMaPzOa5PTOzwpTN3UefBX4MTCZ51uBWSZ/JdWCDSURQV59x05GZFb1sagoXAQsiYgeApH8H6oD/yGVgg0n9ph1s2t7qpiMzK3rZ3H1Uz+7JowJ4PjfhDE67+hNqZngQPDMrbtnUFHYAyyTdTTJ20enAfZL+EyAiLsthfINC3aoM1UMrmHPwyHyHYmaWU9kkhf9LX7sszlEsg1bdygwnzBhLWVnJPZ5hZiUmm1tSrx+IQAarxp1tPLO+ibOOmZzvUMzMci6bPoWS9ujqLURAzUx3MptZ8XNS2IO6+gxlguOmeU5mMyt+WScFSSX51FZd/WaOmDSKkUOzfc7PzKxwZfPw2kJJjwPPpuvHSfqvnEc2CLR3dPLoqi1+PsHMSkY2NYUrSeZn3gQQEY8Bb8hlUIPF8pea2N7a4f4EMysZ2SSFsoio77GtIxfBDDZdM61Nd1Iws9KQTUP5akkLgZBUDnwceCa3YQ0OtfUZDq4eytSxw/IdipnZgMimpvAR4DKSqThfAk5KtxW9uvoMC2aMRfJDa2ZWGrJ5eG09cMEAxDKovLS1mYbMTi4+eWa+QzEzGzB7TAqSriUZ82g3EXFJTiIaJJZ4Uh0zK0HZ9Cn8rttyFfB2YHVuwhk8auszVFaUcdQho/MdipnZgMmm+ej27uuSbgF+m7OIBom6+gzHTR1NZYUf+jaz0rEvn3iHAjMOdCCDSXNbB8tebGSB508wsxKTTZ9Chpf7FMqAzcDluQwq3x5f00hbR7g/wcxKTr9JQcm9mMcBa9JNnRHxik7nYlO7ctdDax4Ez8xKS7/NR2kC+N+I6EhfRZ8QIOlPOHT8CA4aWZJjAJpZCcumT+EhSfNzHskgEREsWZVx05GZlaQ+m48kVUREO/Aa4MOSngO2AyKpRBRloli5aQebt7c6KZhZSeqvT+EhYD7wtgGKZVCoXbkZ8ENrZlaa+ksKAoiI5wYolkFhyaoMo6oqmD1hZL5DMTMbcP0lhQmSLutrZ0T8554OLukM4LtAOXBdRHy1x/7LgA8B7cAG4AO9DNM9oOrqM8yfMZayMg+CZ2alp7+O5nJgJFDdx6tf6TDbVwFnAvOACyXN61HsEaAmIo4Ffgp8fW8v4EBq3NnGMy9tY4HnTzCzEtVfTWFtRFyxH8deCKyIiOcBJN0GnAs8uatARNzbrfxi4KL9ON9+W7LKg+CZWWnrr6awv+0nU9h94LyGdFtfPgj8qtdApEsk1Uqq3bBhw36G1bcl9RnKy8Rx0/zQmpmVpv6Swmn7eezekkqvD79JugioAb7R2/6IuCYiaiKiZsKECfsZVt/q6jMcObmaEUOzGTzWzKz49JkUImLzfh67AZjWbX0q8GLPQpLeBHwWOCciWvbznPusvaOTR1dvcX+CmZW0XI4L/TAwR9KhkipJZm9b1L2ApBOAH5IkhPU5jGWPnl7XxI7WDua7P8HMSljOkkL6NPTHgLuBp4A7ImKZpCsknZMW+wbJHU4/kfSopEV9HC7n6jzTmplZVjOv7bOIuAu4q8e2L3RbflMuz7836uozTBpVxZQxw/IdiplZ3nhasVRdfTIIXjJauJlZaXJSANY1NrNmy073J5hZyXNSwP0JZma7OCmQJIWqIWUcdciofIdiZpZXTgpA3aoMx04dw5By/3OYWWkr+U/Bna0dLFvT6KYjMzOcFFjasIX2zvCTzGZmOClQl46M6juPzMycFFhSn+GwCSMYN6Iy36GYmeVdSSfJAn0jAAALl0lEQVSFiEgeWnPTkZkZUOJJ4fmN28nsaHMns5lZqqSTwq6H1mpmOimYmUGJJ4Ul9RlGDxvCYeNH5jsUM7NBoaSTQm19hvnTx1BW5kHwzMyghJPClh2trFi/zf0JZmbdlGxSeGTVFsDPJ5iZdVeySaGuPkN5mTh+2ph8h2JmNmiUbFKord/MvMmjGF6Z08nnzMwKSkkmhbaOTh5b7UHwzMx6Ksmk8PTaJna2dbg/wcysh5JMCnX1mwGocVIwM9tNSSaF2voMk0dXcciYYfkOxcxsUCnJpLCkPuOmIzOzXpRcUnhxy05ebGz2yKhmZr0ouaSwZJUHwTMz60vJJYXalRmqhpRx5ORR+Q7FzGzQKbmksGRVhuOmjmFIecldupnZHpXUJ+OO1naWvbjVD62ZmfWhpJLC0oZGOjrD/QlmZn0oqaSwa6a1E6Y5KZiZ9abkksKsCSMYO6Iy36GYmQ1KOU0Kks6QtFzSCkmX97J/qKTb0/0PSpqZq1g6O4MlqzLuTzAz60fOkoKkcuAq4ExgHnChpHk9in0QyETEbODbwNdyFc/zG7ezZUcbNTPG5eoUZmYFL5c1hYXAioh4PiJagduAc3uUORe4KV3+KXCapJxMmLxrEDwPb2Fm1rdcJoUpwOpu6w3ptl7LREQ70Agc1PNAki6RVCupdsOGDfsUzNjhlbx53kQOGz9in95vZlYKcjntWG/f+GMfyhAR1wDXANTU1LxifzZOP2oSpx81aV/eamZWMnJZU2gApnVbnwq82FcZSRXAaGBzDmMyM7N+5DIpPAzMkXSopErgAmBRjzKLgPeny+cBv4+IfaoJmJnZ/stZ81FEtEv6GHA3UA7cEBHLJF0B1EbEIuB64BZJK0hqCBfkKh4zM9uzXPYpEBF3AXf12PaFbsvNwDtzGYOZmWWvpJ5oNjOz/jkpmJlZFycFMzPr4qRgZmZdVGh3gEraANTv49vHAxsPYDiFpFSv3dddWnzdfZsRERP2dKCCSwr7Q1JtRNTkO458KNVr93WXFl/3/nPzkZmZdXFSMDOzLqWWFK7JdwB5VKrX7usuLb7u/VRSfQpmZta/UqspmJlZP5wUzMysS1EmBUlnSFouaYWky3vZP1TS7en+ByXNHPgoD7wsrvsySU9KWirpHkkz8hHngban6+5W7jxJIaloblnM5tolvSv9vS+TdOtAx5gLWfytT5d0r6RH0r/3s/IR54Em6QZJ6yU90cd+Sboy/XdZKmn+Xp8kIorqRTJM93PAYUAl8Bgwr0eZvwOuTpcvAG7Pd9wDdN1vAIanyx8pletOy1UDfwIWAzX5jnsAf+dzgEeAsen6wfmOe4Cu+xrgI+nyPGBlvuM+QNf+OmA+8EQf+88CfkUyq+VJwIN7e45irCksBFZExPMR0QrcBpzbo8y5wE3p8k+B0yT1NjVoIdnjdUfEvRGxI11dTDIbXqHL5vcN8CXg60DzQAaXY9lc+4eBqyIiAxAR6wc4xlzI5roDGJUuj+aVsz4WpIj4E/3PTnkucHMkFgNjJE3em3MUY1KYAqzutt6Qbuu1TES0A43AQQMSXe5kc93dfZDkG0Wh2+N1SzoBmBYRdw5kYAMgm9/5XGCupPslLZZ0xoBFlzvZXPcXgYskNZDM6fLxgQkt7/b2c+AVcjrJTp709o2/53232ZQpNFlfk6SLgBrg9TmNaGD0e92SyoBvAxcPVEADKJvfeQVJE9KpJDXDP0s6OiK25Di2XMrmui8EboyIb0l6NckMj0dHRGfuw8ur/f5sK8aaQgMwrdv6VF5ZdewqI6mCpHrZX5WsEGRz3Uh6E/BZ4JyIaBmg2HJpT9ddDRwN/EHSSpJ21kVF0tmc7d/6LyKiLSJeAJaTJIlCls11fxC4AyAiHgCqSAaNK3ZZfQ70pxiTwsPAHEmHSqok6Uhe1KPMIuD96fJ5wO8j7aUpYHu87rQZ5YckCaEY2pZhD9cdEY0RMT4iZkbETJK+lHMiojY/4R5Q2fyt/5zkBgMkjSdpTnp+QKM88LK57lXAaQCSjiRJChsGNMr8WAS8L70L6SSgMSLW7s0Biq75KCLaJX0MuJvkLoUbImKZpCuA2ohYBFxPUp1cQVJDuCB/ER8YWV73N4CRwE/SfvVVEXFO3oI+ALK87qKU5bXfDZwu6UmgA/jHiNiUv6j3X5bX/WngWkmfImk+ubgIvvgh6cckTYHj0/6SfwWGAETE1ST9J2cBK4AdwN/s9TmK4N/JzMwOkGJsPjIzs33kpGBmZl2cFMzMrIuTgpmZdXFSMDOzLk4KllOSOiQ92u01s5+yM/sa/XGgSaqRdGW6fKqkk7vtu1TS+wYwluP3ZZRPSZMl3Zkunyqpsdvv4Xfp9i9KWpNue0LSOb1sf1LShd2O+01JbzxQ12eDS9E9p2CDzs6IOD7fQeyt9OG2XQ+4nQpsA/6S7rv6QJ9PUkU6DldvjicZluSuvTzsZcC13db/HBF/1Uu5b0fEN9OHvP4s6eAe2+cAdZJ+GhFtwH+lx/39XsZjBcA1BRtwaY3gz5KWpK+TeylzlKSH0m+qS9MPJiRd1G37DyWV9/LelZK+lpZ7SNLsdPsMJfNI7JpPYnq6/Z3pt+THJP0p3XaqpDvTms2lwKfSc742/Rb9D5KOlPRQj+tami4vkPRHSXWS7lYvI1VKulHSf0q6F/iapIWS/qJkDoC/SDo8fWL3CuD89PznSxqhZFz9h9OyvY0KC/DXwK+z/b1ExFNAOz2Gg4iIZ0kehBqbrtcDB0malO2xrXA4KViuDevWZPG/6bb1wJsjYj5wPnBlL++7FPhuWsuoARrSb7LnA6ek2zuA9/Rx3q0RsRD4HvCddNv3SIYVPhb4727n/QLwlog4DtjtCe+IWAlcTfKt+fiI+HO3fU8BlZIOSzedD9whaQjJt+nzImIBcAPw733EORd4U0R8GngaeF1EnJDG9JV0aOgvkMx9cXxE3E4ydtXvI+JEkiEsviFpRPeDSjoUyPQY3+q13X4Xn+0ZiKRXAZ30GA5CyUQtz/YYGmUJcEof12QFzM1Hlmu9NR8NAb4nadcH+9xe3vcA8FlJU4GfRcSzkk4DFgAPKxmmYxhJgunNj7v9/Ha6/GrgHenyLSTzKwDcD9wo6Q7gZ3tzcSSDrr0L+CpJUjgfOJxkEL7fpnGWA32NP/OTiOhIl0cDN6W1oiAdvqAXpwPnSPqHdL0KmA481a3MZF451k9fzUefUjJybhNwfkREGvenJH2YZDKbnkNurwcO6SM+K2BOCpYPnwJeAo4jqa2+YuKbiLhV0oPA2cDdkj5EMizwTRHxmSzOEX0sv6JMRFyafks+G3g0TVbZup1kLKmfJYeKZyUdAyyLiFdn8f7t3Za/BNwbEW9Pm63+0Md7BPx1RCzv57g7SZJFNr4dEd/sa7ukdwA3S5oVEbt+V1XpOazIuPnI8mE0sDYd2/69JN+kd5M2yTwfEVeSjPx4LHAPcN6ujlBJ49T3PNPnd/v5QLr8F14e/PA9wH3pcWZFxIMR8QVgI7sPPQzJN+jq3k4SEc+R1HY+T5IgIBmeeoKScfyRNETSUX3E2d1oYE26fHE/578b+LjSr/NKRr/t6RlgZhbn3KOI+BlJp/v7u22eCwyKO8XswHJSsHz4PvB+SYtJPly291LmfOAJSY8CR5D0BTwJfA74Tdqh+1uSZpLeDE1rGn9PUjMB+ATwN+l735vug6RN/nElt8P+iWTO3+5+Cbx9V0dzL+e6HbiIl8fvbyUZkv1rkh4DHgVe0Znei68D/yHpfnZPlPcC83Z1NJPUKIYAS9OYv9TzQBGxHXhuVyf7AXAFcJmksrTPZDYv351lRcSjpFrRUTKZTk1EbMx3LPkk6e3Agoj4XA6OOz8iPn8gj2uDg/sUzIpURPyvpFzMPV4BfCsHx7VBwDUFMzPr4j4FMzPr4qRgZmZdnBTMzKyLk4KZmXVxUjAzsy7/H5YBz0IrFkGlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(rf[\"FPR\"], rf[\"TPR\"])\n",
    "ax = plt.gca()\n",
    "ax.set_xlabel(\"False positive rate (FPR)\")\n",
    "ax.set_ylabel(\"True positive rate (TPR)\")\n",
    "ax.set_title(\"ROC curve\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Threshold</th>\n",
       "      <th>TP</th>\n",
       "      <th>TN</th>\n",
       "      <th>FP</th>\n",
       "      <th>FN</th>\n",
       "      <th>FPR</th>\n",
       "      <th>TPR</th>\n",
       "      <th>ACC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2182</td>\n",
       "      <td>0</td>\n",
       "      <td>1090</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.666870</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>100000</td>\n",
       "      <td>691</td>\n",
       "      <td>1982</td>\n",
       "      <td>200</td>\n",
       "      <td>399</td>\n",
       "      <td>0.091659</td>\n",
       "      <td>0.633945</td>\n",
       "      <td>0.816932</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>200000</td>\n",
       "      <td>859</td>\n",
       "      <td>1764</td>\n",
       "      <td>418</td>\n",
       "      <td>231</td>\n",
       "      <td>0.191567</td>\n",
       "      <td>0.788073</td>\n",
       "      <td>0.801650</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>300000</td>\n",
       "      <td>952</td>\n",
       "      <td>1556</td>\n",
       "      <td>626</td>\n",
       "      <td>138</td>\n",
       "      <td>0.286893</td>\n",
       "      <td>0.873394</td>\n",
       "      <td>0.766504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>400000</td>\n",
       "      <td>997</td>\n",
       "      <td>1339</td>\n",
       "      <td>843</td>\n",
       "      <td>93</td>\n",
       "      <td>0.386343</td>\n",
       "      <td>0.914679</td>\n",
       "      <td>0.713936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>500000</td>\n",
       "      <td>1029</td>\n",
       "      <td>1162</td>\n",
       "      <td>1020</td>\n",
       "      <td>61</td>\n",
       "      <td>0.467461</td>\n",
       "      <td>0.944037</td>\n",
       "      <td>0.669621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>600000</td>\n",
       "      <td>1059</td>\n",
       "      <td>976</td>\n",
       "      <td>1206</td>\n",
       "      <td>31</td>\n",
       "      <td>0.552704</td>\n",
       "      <td>0.971560</td>\n",
       "      <td>0.621944</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>700000</td>\n",
       "      <td>1073</td>\n",
       "      <td>872</td>\n",
       "      <td>1310</td>\n",
       "      <td>17</td>\n",
       "      <td>0.600367</td>\n",
       "      <td>0.984404</td>\n",
       "      <td>0.594438</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>800000</td>\n",
       "      <td>1076</td>\n",
       "      <td>762</td>\n",
       "      <td>1420</td>\n",
       "      <td>14</td>\n",
       "      <td>0.650779</td>\n",
       "      <td>0.987156</td>\n",
       "      <td>0.561736</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>900000</td>\n",
       "      <td>1080</td>\n",
       "      <td>655</td>\n",
       "      <td>1527</td>\n",
       "      <td>10</td>\n",
       "      <td>0.699817</td>\n",
       "      <td>0.990826</td>\n",
       "      <td>0.530257</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1000000</td>\n",
       "      <td>1082</td>\n",
       "      <td>594</td>\n",
       "      <td>1588</td>\n",
       "      <td>8</td>\n",
       "      <td>0.727773</td>\n",
       "      <td>0.992661</td>\n",
       "      <td>0.512225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>1100000</td>\n",
       "      <td>1086</td>\n",
       "      <td>498</td>\n",
       "      <td>1684</td>\n",
       "      <td>4</td>\n",
       "      <td>0.771769</td>\n",
       "      <td>0.996330</td>\n",
       "      <td>0.484108</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1200000</td>\n",
       "      <td>1088</td>\n",
       "      <td>429</td>\n",
       "      <td>1753</td>\n",
       "      <td>2</td>\n",
       "      <td>0.803391</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.463631</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1300000</td>\n",
       "      <td>1088</td>\n",
       "      <td>399</td>\n",
       "      <td>1783</td>\n",
       "      <td>2</td>\n",
       "      <td>0.817140</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.454462</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1400000</td>\n",
       "      <td>1088</td>\n",
       "      <td>351</td>\n",
       "      <td>1831</td>\n",
       "      <td>2</td>\n",
       "      <td>0.839138</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.439792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1500000</td>\n",
       "      <td>1088</td>\n",
       "      <td>292</td>\n",
       "      <td>1890</td>\n",
       "      <td>2</td>\n",
       "      <td>0.866178</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.421760</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>1600000</td>\n",
       "      <td>1088</td>\n",
       "      <td>266</td>\n",
       "      <td>1916</td>\n",
       "      <td>2</td>\n",
       "      <td>0.878093</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.413814</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>1700000</td>\n",
       "      <td>1088</td>\n",
       "      <td>220</td>\n",
       "      <td>1962</td>\n",
       "      <td>2</td>\n",
       "      <td>0.899175</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.399756</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>1800000</td>\n",
       "      <td>1088</td>\n",
       "      <td>197</td>\n",
       "      <td>1985</td>\n",
       "      <td>2</td>\n",
       "      <td>0.909716</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.392726</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>1900000</td>\n",
       "      <td>1088</td>\n",
       "      <td>179</td>\n",
       "      <td>2003</td>\n",
       "      <td>2</td>\n",
       "      <td>0.917965</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.387225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>2000000</td>\n",
       "      <td>1088</td>\n",
       "      <td>160</td>\n",
       "      <td>2022</td>\n",
       "      <td>2</td>\n",
       "      <td>0.926673</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.381418</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>2100000</td>\n",
       "      <td>1088</td>\n",
       "      <td>130</td>\n",
       "      <td>2052</td>\n",
       "      <td>2</td>\n",
       "      <td>0.940422</td>\n",
       "      <td>0.998165</td>\n",
       "      <td>0.372249</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>2200000</td>\n",
       "      <td>1089</td>\n",
       "      <td>110</td>\n",
       "      <td>2072</td>\n",
       "      <td>1</td>\n",
       "      <td>0.949588</td>\n",
       "      <td>0.999083</td>\n",
       "      <td>0.366443</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>2300000</td>\n",
       "      <td>1089</td>\n",
       "      <td>100</td>\n",
       "      <td>2082</td>\n",
       "      <td>1</td>\n",
       "      <td>0.954170</td>\n",
       "      <td>0.999083</td>\n",
       "      <td>0.363386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>2400000</td>\n",
       "      <td>1089</td>\n",
       "      <td>86</td>\n",
       "      <td>2096</td>\n",
       "      <td>1</td>\n",
       "      <td>0.960587</td>\n",
       "      <td>0.999083</td>\n",
       "      <td>0.359108</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>2500000</td>\n",
       "      <td>1089</td>\n",
       "      <td>73</td>\n",
       "      <td>2109</td>\n",
       "      <td>1</td>\n",
       "      <td>0.966544</td>\n",
       "      <td>0.999083</td>\n",
       "      <td>0.355134</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>2600000</td>\n",
       "      <td>1089</td>\n",
       "      <td>73</td>\n",
       "      <td>2109</td>\n",
       "      <td>1</td>\n",
       "      <td>0.966544</td>\n",
       "      <td>0.999083</td>\n",
       "      <td>0.355134</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>2700000</td>\n",
       "      <td>1090</td>\n",
       "      <td>59</td>\n",
       "      <td>2123</td>\n",
       "      <td>0</td>\n",
       "      <td>0.972961</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.351161</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>2800000</td>\n",
       "      <td>1090</td>\n",
       "      <td>55</td>\n",
       "      <td>2127</td>\n",
       "      <td>0</td>\n",
       "      <td>0.974794</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.349939</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>2900000</td>\n",
       "      <td>1090</td>\n",
       "      <td>47</td>\n",
       "      <td>2135</td>\n",
       "      <td>0</td>\n",
       "      <td>0.978460</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.347494</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Threshold    TP    TN    FP    FN       FPR       TPR       ACC\n",
       "0           0     0  2182     0  1090  0.000000  0.000000  0.666870\n",
       "1      100000   691  1982   200   399  0.091659  0.633945  0.816932\n",
       "2      200000   859  1764   418   231  0.191567  0.788073  0.801650\n",
       "3      300000   952  1556   626   138  0.286893  0.873394  0.766504\n",
       "4      400000   997  1339   843    93  0.386343  0.914679  0.713936\n",
       "5      500000  1029  1162  1020    61  0.467461  0.944037  0.669621\n",
       "6      600000  1059   976  1206    31  0.552704  0.971560  0.621944\n",
       "7      700000  1073   872  1310    17  0.600367  0.984404  0.594438\n",
       "8      800000  1076   762  1420    14  0.650779  0.987156  0.561736\n",
       "9      900000  1080   655  1527    10  0.699817  0.990826  0.530257\n",
       "10    1000000  1082   594  1588     8  0.727773  0.992661  0.512225\n",
       "11    1100000  1086   498  1684     4  0.771769  0.996330  0.484108\n",
       "12    1200000  1088   429  1753     2  0.803391  0.998165  0.463631\n",
       "13    1300000  1088   399  1783     2  0.817140  0.998165  0.454462\n",
       "14    1400000  1088   351  1831     2  0.839138  0.998165  0.439792\n",
       "15    1500000  1088   292  1890     2  0.866178  0.998165  0.421760\n",
       "16    1600000  1088   266  1916     2  0.878093  0.998165  0.413814\n",
       "17    1700000  1088   220  1962     2  0.899175  0.998165  0.399756\n",
       "18    1800000  1088   197  1985     2  0.909716  0.998165  0.392726\n",
       "19    1900000  1088   179  2003     2  0.917965  0.998165  0.387225\n",
       "20    2000000  1088   160  2022     2  0.926673  0.998165  0.381418\n",
       "21    2100000  1088   130  2052     2  0.940422  0.998165  0.372249\n",
       "22    2200000  1089   110  2072     1  0.949588  0.999083  0.366443\n",
       "23    2300000  1089   100  2082     1  0.954170  0.999083  0.363386\n",
       "24    2400000  1089    86  2096     1  0.960587  0.999083  0.359108\n",
       "25    2500000  1089    73  2109     1  0.966544  0.999083  0.355134\n",
       "26    2600000  1089    73  2109     1  0.966544  0.999083  0.355134\n",
       "27    2700000  1090    59  2123     0  0.972961  1.000000  0.351161\n",
       "28    2800000  1090    55  2127     0  0.974794  1.000000  0.349939\n",
       "29    2900000  1090    47  2135     0  0.978460  1.000000  0.347494"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Customized edge selector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several issues with just thresholding `MLEDistance` (say `MLEDistance < 512Kb`):\n",
    "- Some nodes may get excluded since all their edges have a predicted long distance\n",
    "- Some nodes, particularly contigs with low link density, tend to have small `MLEDistance` but in fact are not at all significant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ground truth path\n",
    "\n",
    "Make a ground truth path and use `allhic.anchor()` to generate matrix and plot it in `jcvi.assembly.hic.heatmap()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "idcChr1\t1\t127964\tidcChr1.ctg1\n",
      "idcChr1\t127965\t176555\tidcChr1.ctg2\n",
      "idcChr1\t176556\t296400\tidcChr1.ctg3\n",
      "idcChr1\t296401\t481842\tidcChr1.ctg4\n",
      "idcChr1\t481843\t725279\tidcChr1.ctg5\n",
      "idcChr1\t725280\t726201\tidcChr1.ctg6\n",
      "idcChr1\t726202\t811858\tidcChr1.ctg7\n",
      "idcChr1\t811859\t1016514\tidcChr1.ctg8\n",
      "idcChr1\t1016515\t1090313\tidcChr1.ctg9\n",
      "idcChr1\t1090314\t1122899\tidcChr1.ctg10\n"
     ]
    }
   ],
   "source": [
    "!head ctg.posi.bed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[0;33m19:37:11 [base]\u001b[0m\u001b[0;35m Load file `ctg.posi.bed`\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "from jcvi.formats.bed import Bed\n",
    "\n",
    "bed = Bed(\"ctg.posi.bed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "contigs = []\n",
    "sizes = []\n",
    "chrs = [\"idcChr1\", \"jpcChr1\"]\n",
    "for c in chrs:\n",
    "    for b in bed.sub_bed(c):\n",
    "        contigs.append(b.accn)\n",
    "        sizes.append(b.span)\n",
    "\n",
    "with open(\"rice_truth.tour\", \"w\") as fw:\n",
    "    print >> fw, \">TRUTH\"\n",
    "    print >> fw, \" \".join(x + \"+\" for x in contigs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "992\n"
     ]
    }
   ],
   "source": [
    "print(len(contigs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([  358346,   709598,  1058981,  1379910,  1682293,  1984208,\n",
      "        2277372,  2567174,  2854846,  3137570,  3412449,  3687014,\n",
      "        3959758,  4232218,  4504105,  4775230,  5046122,  5316991,\n",
      "        5586095,  5854596,  6121265,  6387083,  6649994,  6911770,\n",
      "        7172076,  7431252,  7688990,  7946649,  8201200,  8449286,\n",
      "        8695733,  8939169,  9181289,  9422272,  9662279,  9900212,\n",
      "       10136081, 10370442, 10604155, 10836172, 11067780, 11295635,\n",
      "       11522776, 11745981, 11967507, 12188544, 12408901, 12629167,\n",
      "       12849282, 13069143, 13288070, 13506362, 13724171, 13941773,\n",
      "       14157706, 14373497, 14588504, 14803374, 15016348, 15228971,\n",
      "       15441255, 15652783, 15862866, 16072070, 16278708, 16484828,\n",
      "       16689483, 16893724, 17097495, 17300586, 17503523, 17706253,\n",
      "       17908045, 18109386, 18310683, 18510644, 18709097, 18905683,\n",
      "       19101759, 19297727, 19493625, 19688957, 19883543, 20077142,\n",
      "       20270328, 20463403, 20655499, 20846344, 21035293, 21222585,\n",
      "       21409376, 21596123, 21782287, 21968081, 22153522, 22338842,\n",
      "       22523638, 22708323, 22891562, 23074550, 23257389, 23439209,\n",
      "       23619590, 23798439, 23977139, 24155720, 24333902, 24511931,\n",
      "       24689874, 24867813, 25045594, 25222436, 25398861, 25575040,\n",
      "       25751205, 25926597, 26100866, 26273617, 26445758, 26617017,\n",
      "       26788027, 26958861, 27129383, 27299773, 27469005, 27638200,\n",
      "       27806642, 27974909, 28142810, 28310294, 28477458, 28644152,\n",
      "       28810791, 28977348, 29143722, 29309932, 29475889, 29641504,\n",
      "       29806560, 29969673, 30131622, 30292716, 30453481, 30613897,\n",
      "       30774010, 30934025, 31093484, 31252075, 31410488, 31568595,\n",
      "       31726578, 31884530, 32042258, 32198723, 32355112, 32511167,\n",
      "       32666749, 32822162, 32977532, 33132866, 33287735, 33442063,\n",
      "       33596255, 33749863, 33903261, 34056518, 34209359, 34362155,\n",
      "       34514902, 34666917, 34818807, 34970185, 35120855, 35271424,\n",
      "       35421973, 35572294, 35722443, 35872409, 36021899, 36171353,\n",
      "       36320732, 36470067, 36619267, 36768212, 36917114, 37065755,\n",
      "       37214145, 37362054, 37509646, 37657159, 37804480, 37951456,\n",
      "       38098315, 38244890, 38391122, 38537256, 38683321, 38829084,\n",
      "       38974719, 39120246, 39265479, 39410703, 39555926, 39701042,\n",
      "       39846100, 39990378, 40134584, 40278623, 40422395, 40566001,\n",
      "       40709595, 40852700, 40995599, 41138219, 41280816, 41423293,\n",
      "       41565657, 41707796, 41848745, 41989142, 42129501, 42269766,\n",
      "       42410009, 42549884, 42689629, 42829225, 42968448, 43107545,\n",
      "       43246607, 43384796, 43522137, 43659163, 43795666, 43931782,\n",
      "       44067387, 44202879, 44337888, 44472384, 44606578, 44740674,\n",
      "       44874493, 45008173, 45141669, 45275036, 45407168, 45539114,\n",
      "       45670819, 45802469, 45933645, 46064328, 46193857, 46323244,\n",
      "       46452043, 46580793, 46709403, 46837583, 46965621, 47093584,\n",
      "       47221217, 47348380, 47474613, 47600292, 47725560, 47850510,\n",
      "       47975283, 48099850, 48224157, 48348329, 48472315, 48596012,\n",
      "       48719391, 48842322, 48965043, 49087373, 49209540, 49330922,\n",
      "       49452188, 49573432, 49694257, 49815034, 49935491, 50055640,\n",
      "       50175484, 50295108, 50414507, 50533392, 50652272, 50770905,\n",
      "       50889467, 51008017, 51126278, 51244171, 51361909, 51479458,\n",
      "       51596978, 51714184, 51831073, 51947810, 52064268, 52180706,\n",
      "       52296012, 52411179, 52525861, 52640478, 52755017, 52869292,\n",
      "       52983533, 53097659, 53211384, 53324905, 53438145, 53551369,\n",
      "       53664565, 53777645, 53890398, 54002889, 54115256, 54227109,\n",
      "       54338812, 54449917, 54561008, 54671648, 54781708, 54891694,\n",
      "       55001168, 55110076, 55218881, 55327360, 55435791, 55543930,\n",
      "       55651817, 55759369, 55866855, 55974202, 56081543, 56188651,\n",
      "       56295433, 56401867, 56508227, 56614403, 56720457, 56826137,\n",
      "       56931529, 57036803, 57141894, 57246732, 57351441, 57455944,\n",
      "       57560279, 57664602, 57768923, 57873221, 57977496, 58081754,\n",
      "       58185877, 58289784, 58393154, 58496468, 58599311, 58702119,\n",
      "       58804764, 58907396, 59009135, 59109421, 59209286, 59308365,\n",
      "       59406436, 59504351, 59601850, 59699242, 59796244, 59893092,\n",
      "       59989713, 60086255, 60182518, 60278520, 60374499, 60470073,\n",
      "       60565401, 60660511, 60755562, 60850074, 60944508, 61038841,\n",
      "       61132885, 61226902, 61320746, 61414530, 61508284, 61601438,\n",
      "       61694063, 61786332, 61878335, 61970306, 62062274, 62154208,\n",
      "       62245998, 62336993, 62427689, 62518293, 62608875, 62699294,\n",
      "       62789645, 62879988, 62970318, 63060564, 63150704, 63240436,\n",
      "       63329959, 63419452, 63508877, 63598220, 63687424, 63776450,\n",
      "       63865374, 63954233, 64042948, 64131641, 64220048, 64308441,\n",
      "       64396833, 64484921, 64572899, 64660711, 64748317, 64835921,\n",
      "       64923369, 65010765, 65097962, 65185096, 65271793, 65358473,\n",
      "       65444804, 65531044, 65616752, 65702408, 65787786, 65873134,\n",
      "       65958426, 66043640, 66128595, 66213147, 66297272, 66380962,\n",
      "       66464522, 66548050, 66631378, 66713985, 66796500, 66878977,\n",
      "       66961252, 67043354, 67125333, 67206719, 67288022, 67369242,\n",
      "       67450320, 67531103, 67611804, 67691912, 67771883, 67851669,\n",
      "       67931285, 68010676, 68089869, 68168938, 68247836, 68326501,\n",
      "       68405127, 68483308, 68561336, 68639253, 68716867, 68794061,\n",
      "       68871155, 68948144, 69025133, 69101844, 69178536, 69254917,\n",
      "       69331237, 69407526, 69483631, 69559721, 69635703, 69711450,\n",
      "       69786757, 69861396, 69935949, 70010374, 70084669, 70158915,\n",
      "       70232713, 70306233, 70379653, 70452927, 70525898, 70598861,\n",
      "       70671812, 70744687, 70817444, 70889903, 70962336, 71034718,\n",
      "       71107031, 71179131, 71251142, 71323076, 71394981, 71466681,\n",
      "       71537927, 71609165, 71680363, 71751461, 71822332, 71893136,\n",
      "       71963827, 72034368, 72104880, 72175384, 72245608, 72315791,\n",
      "       72385966, 72455602, 72524933, 72593610, 72662225, 72730695,\n",
      "       72798783, 72866865, 72934827, 73002676, 73070472, 73138152,\n",
      "       73205720, 73272800, 73339860, 73406886, 73473841, 73540552,\n",
      "       73607171, 73673511, 73739542, 73805550, 73871536, 73937336,\n",
      "       74003091, 74068625, 74134121, 74199559, 74264921, 74330028,\n",
      "       74394693, 74459219, 74523571, 74587727, 74651433, 74715116,\n",
      "       74778792, 74842352, 74905848, 74969297, 75032695, 75095907,\n",
      "       75159100, 75222245, 75285016, 75347763, 75410393, 75472875,\n",
      "       75535196, 75597345, 75659440, 75721453, 75783137, 75844600,\n",
      "       75906000, 75967180, 76028133, 76089030, 76149899, 76210760,\n",
      "       76271324, 76331794, 76392102, 76452387, 76512423, 76572333,\n",
      "       76632111, 76691724, 76751320, 76810857, 76870250, 76929482,\n",
      "       76988711, 77047783, 77106211, 77164561, 77222873, 77280925,\n",
      "       77338947, 77396757, 77454543, 77512295, 77569905, 77627236,\n",
      "       77684313, 77741259, 77797907, 77854432, 77910777, 77967018,\n",
      "       78023173, 78079315, 78135222, 78191101, 78246795, 78302435,\n",
      "       78358069, 78413575, 78468920, 78524155, 78579303, 78634278,\n",
      "       78689253, 78743817, 78798232, 78852444, 78906572, 78960134,\n",
      "       79013664, 79067061, 79120411, 79173655, 79226692, 79279422,\n",
      "       79331841, 79384103, 79436205, 79488117, 79539984, 79591551,\n",
      "       79643050, 79694041, 79744702, 79795277, 79845832, 79896140,\n",
      "       79946377, 79996608, 80046805, 80096974, 80147005, 80196685,\n",
      "       80246093, 80295431, 80344666, 80393892, 80442624, 80491244,\n",
      "       80539834, 80588423, 80636898, 80685080, 80733222, 80781070,\n",
      "       80828492, 80875809, 80923100, 80970205, 81017253, 81064050,\n",
      "       81110800, 81157508, 81204105, 81250500, 81296779, 81343017,\n",
      "       81389156, 81435160, 81480942, 81526590, 81572043, 81617367,\n",
      "       81662501, 81707626, 81752593, 81797351, 81841943, 81886441,\n",
      "       81930674, 81974810, 82018912, 82062923, 82106514, 82149786,\n",
      "       82192677, 82235371, 82277617, 82319814, 82361929, 82404033,\n",
      "       82445883, 82487487, 82529084, 82570673, 82612187, 82653676,\n",
      "       82695053, 82736071, 82776979, 82817873, 82858765, 82899309,\n",
      "       82939822, 82980199, 83020459, 83060629, 83100782, 83140759,\n",
      "       83180642, 83220480, 83260315, 83300097, 83339845, 83379224,\n",
      "       83418193, 83457085, 83495977, 83534456, 83572552, 83610564,\n",
      "       83648163, 83685623, 83723072, 83760391, 83797701, 83834919,\n",
      "       83872109, 83909272, 83946413, 83983187, 84019961, 84056660,\n",
      "       84093350, 84129989, 84166186, 84202245, 84238279, 84274257,\n",
      "       84310076, 84345401, 84380646, 84415769, 84450750, 84485672,\n",
      "       84520222, 84554698, 84589106, 84623168, 84657178, 84691135,\n",
      "       84724856, 84758542, 84792019, 84825049, 84857898, 84890680,\n",
      "       84923265, 84955666, 84987929, 85020171, 85052256, 85084123,\n",
      "       85115536, 85146731, 85177746, 85208757, 85239706, 85270344,\n",
      "       85300965, 85331492, 85361986, 85392437, 85422805, 85453059,\n",
      "       85483292, 85513143, 85542724, 85572064, 85601396, 85630665,\n",
      "       85659885, 85689075, 85718142, 85747196, 85776224, 85804965,\n",
      "       85833513, 85861985, 85890047, 85918037, 85945920, 85973525,\n",
      "       86001048, 86028105, 86054983, 86081788, 86108525, 86135158,\n",
      "       86161769, 86187988, 86214120, 86239833, 86265322, 86290527,\n",
      "       86315719, 86340821, 86365687, 86390353, 86415014, 86439604,\n",
      "       86464182, 86488638, 86513069, 86537433, 86561559, 86585668,\n",
      "       86609581, 86633232, 86656555, 86679532, 86702454, 86725206,\n",
      "       86747919, 86770586, 86793249, 86815875, 86838289, 86860682,\n",
      "       86883050, 86905282, 86927398, 86949206, 86970994, 86992629,\n",
      "       87014225, 87035784, 87057220, 87078559, 87099803, 87120858,\n",
      "       87141713, 87162075, 87182424, 87202667, 87222830, 87242693,\n",
      "       87262429, 87282114, 87301750, 87321290, 87340637, 87359570,\n",
      "       87378366, 87396970, 87415569, 87433968, 87452198, 87470421,\n",
      "       87488523, 87506578, 87524628, 87542592, 87560554, 87578463,\n",
      "       87596298, 87614099, 87631700, 87649284, 87666867, 87684393,\n",
      "       87701900, 87719069, 87736116, 87753141, 87770143, 87786653,\n",
      "       87803098, 87819456, 87835718, 87851742, 87867551, 87883314,\n",
      "       87899055, 87914667, 87930160, 87945561, 87960953, 87976037,\n",
      "       87990920, 88005646, 88020221, 88034769, 88049269, 88063732,\n",
      "       88077851, 88091536, 88105116, 88118693, 88132023, 88145301,\n",
      "       88158505, 88171453, 88184376, 88197173, 88209893, 88222451,\n",
      "       88234631, 88246588, 88258493, 88270374, 88282130, 88293763,\n",
      "       88304808, 88315837, 88326861, 88337713, 88348522, 88358953,\n",
      "       88369368, 88379609, 88389537, 88399290, 88408971, 88418365,\n",
      "       88427661, 88436898, 88445537, 88454040, 88462237, 88470352,\n",
      "       88478201, 88486033, 88493782, 88501493, 88509185, 88516844,\n",
      "       88524369, 88531763, 88539036, 88546302, 88553131, 88559805,\n",
      "       88566055, 88572041, 88577587, 88582909, 88587994, 88592972,\n",
      "       88597619, 88602195, 88606714, 88610999, 88614985, 88618867,\n",
      "       88622714, 88626540, 88629888, 88633108, 88636164, 88639176,\n",
      "       88642109, 88645023, 88647855, 88650433, 88652966, 88655430,\n",
      "       88657852, 88659975, 88661739, 88663114, 88664286, 88665294,\n",
      "       88666215, 88666873, 88667530, 88668164, 88668789, 88669395,\n",
      "       88669820, 88670167]), 135009, 237)\n"
     ]
    }
   ],
   "source": [
    "from jcvi.assembly.base import calculate_A50\n",
    "\n",
    "print(calculate_A50(sizes))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
